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 26, 2020, 07:22
Default
  #41
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, 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).
But are you aware about the meaning of the FFT? It is not a solver for computing the 5-point FD-based linear Poisson equation.
I think you have not understood that.
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 07:24
Default
  #42
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
But are you aware about the meaning of the FFT? It is not a solver for computing the 5-point FD-based linear Poisson equation.
I think you have not understood that.

Are you sure we're talking about the same FastPoissonSolver ? One like explained here ? http://citeseerx.ist.psu.edu/viewdoc...=rep1&type=pdf
Eifoehn4 likes this.
XorT is offline   Reply With Quote

Old   June 26, 2020, 07:39
Default
  #43
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
Are you sure we're talking about the same FastPoissonSolver ? One like explained here ? http://citeseerx.ist.psu.edu/viewdoc...=rep1&type=pdf
I will read but the term “fast Poisson” does not mean univocally the spectral method based on FFT. Maybe is something like pseudo-spectral... after reading I can tell you ..
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 07:41
Default
  #44
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 will read but the term “fast Poisson” does not mean univocally the spectral method based on FFT. Maybe is something like pseudo-spectral... after reading I can tell you ..

possibly as i said, i'm no expert on this.
i was only mentioning "FFT based FastPoisson solver" initially and the paper also talks about "solving directly using a fast Poisson solver". if you were interpreting this as something else then that would at least explain all this confusion.
XorT is offline   Reply With Quote

Old   June 26, 2020, 08:25
Default
  #45
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
possibly as i said, i'm no expert on this.
i was only mentioning "FFT based FastPoisson solver" initially and the paper also talks about "solving directly using a fast Poisson solver". if you were interpreting this as something else then that would at least explain all this confusion.



Ok, I read and indeed no pure spectral representation of the solution is implied, it is just a DST method to solve the second order accurate FD equation, like SOR or multigrid or GMRES method. I am not sure it is really computationally cheap...
However, now we are sure that the method is congruent to the exact projection of Chorin and is second order in space.

I suggest to test your code by starting with small density differences and then, if it works, increasing the density gap.
Are you using a full explicit method for convection and diffusion? How do you set the time step for the numerical stability constraint?
Eifoehn4 likes this.
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 08:33
Default
  #46
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Ok, I read and indeed no pure spectral representation of the solution is implied, it is just a DST method to solve the second order accurate FD equation, like SOR or multigrid or GMRES method. I am not sure it is really computationally cheap...
However, now we are sure that the method is congruent to the exact projection of Chorin and is second order in space.

Thanks god we sorted this out, haha

Actually it is very cheap compared to all alternatives. Especially on a GPU. My shader implementation needs 0.6ms on a 1024^2 grid (on a GeForce 1070).



Quote:
Originally Posted by FMDenaro View Post
I suggest to test your code by starting with small density differences and then, if it works, increasing the density gap.
Are you using a full explicit method for convection and diffusion? How do you set the time step for the numerical stability constraint?

I use the MacCormack scheme for advection. No explicit diffusion is added, neither are any forces. See my initial post. Even for very small time-steps/velocities of about 1/100 of a gridcell per timestep will make it blow up eventually. Also smaller density ratios will make it blow-up. It will just take a little longer. As stated, an exception is a density ratio of 1 in which case it works like a charm (but in that case the "modified" part of the pressure solver is not actually in effect).
XorT is offline   Reply With Quote

Old   June 26, 2020, 08:38
Default
  #47
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
Thanks god we sorted this out, haha

Actually it is very cheap compared to all alternatives. Especially on a GPU. My shader implementation needs 0.6ms on a 1024^2 grid (on a GeForce 1070).






I use the MacCormack scheme for advection. No explicit diffusion is added, neither are any forces. See my initial post. Even for very small time-steps/velocities of about 1/100 of a gridcell per timestep will make it blow up eventually. Also smaller density ratios will make it blow-up. It will just take a little longer. As stated, an exception is a density ratio of 1 in which case it works like a charm (but in that case the "modified" part of the pressure solver is not actually in effect).

Well, actually you are solving Euler equations with an initial discontinuity. Without diffusion, using central space discretization, you will get strong oscillations and finally instability. Have a look to some shock capturing method based on the artificial dissipation method
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 08:43
Default
  #48
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
Well, actually you are solving Euler equations with an initial discontinuity. Without diffusion, using central space discretization, you will get strong oscillations and finally instability. Have a look to some shock capturing method based on the artificial dissipation method

What speaks against it is that NOT advecting the fluid will avoid the blow up and it stabilizes (with density ratio of 1000) to a nice solution of divergence free velocities (in which the water is mostly unaffected, and the air flows around it to counter divergence).


here's a picture displaying vertical velocity color coded (green positive, red negative): https://ibb.co/GTCDZ81 (left initial condition, right solution after few iterations)
XorT is offline   Reply With Quote

Old   June 26, 2020, 08:56
Default
  #49
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
What speaks against it is that NOT advecting the fluid will avoid the blow up and it stabilizes (with density ratio of 1000) to a nice solution of divergence free velocities (in which the water is mostly unaffected, and the air flows around it to counter divergence).


here's a picture displaying vertical velocity color coded (green positive, red negative): https://ibb.co/GTCDZ81 (left initial condition, right solution after few iterations)
You do not need accuracy, thus I suggest using first order explicit time integration and first order upwind for the convective terms. That would stabilize the solution owing to the numerical diffusion. Just be aware about the cfl condition.
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 09:06
Default
  #50
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
What speaks against it is that NOT advecting the fluid will avoid the blow up and it stabilizes (with density ratio of 1000) to a nice solution of divergence free velocities (in which the water is mostly unaffected, and the air flows around it to counter divergence).

here's a picture displaying vertical velocity color coded (green positive, red negative): https://ibb.co/GTCDZ81 (left initial condition, right solution after few iterations)

Thanks you for the references. Now it is more clear what you are doing. I would follow the hints made by FMDenaro.

Moreover i think the key problem with high density ratios is this term in your modified equation:

\dots = \nabla \cdot \left[  \left( 1 -  \frac{ \rho_{0} }{ \rho^{n+1}} \right) \nabla \hat{p} \right] + \dots

How do you approximate this at the moment? Did you checked that this is correct?

Edit:

Quote:
Originally Posted by FMDenaro View Post
... it is just a DST method to solve the second order accurate FD equation, like SOR or multigrid or GMRES method. I am not sure it is really computationally cheap...
I think a simple vectorized version of the Jacobi method would also do the job, especially if you don't need accuracy. Only one iteration each time step, this would result in O(N) instead of O(N \text{log}(N)).
FMDenaro likes this.
__________________
Check out my side project:

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

Old   June 26, 2020, 09:35
Default
  #51
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
Thanks you for the references. Now it is more clear what you are doing. I would follow the hints made by FMDenaro.

Moreover i think the key problem with high density ratios is this term in your modified equation:

\dots = \nabla \cdot \left[  \left( 1 -  \frac{ \rho_{0} }{ \rho^{n+1}} \right) \nabla \hat{p} \right] + \dots

How do you approximate this at the moment? Did you checked that this is correct?

Edit:



I think a simple vectorized version of the Jacobi method, only one iteration each step, would do the job, especially if you don't need accuracy. The computional cound would be O(N) instead of O(N log(N)).

you mean how do i approximate \hat{p} ? as suggested in the paper as \hat{p}^{n+1}=2p^{n}-p^{n-1}


I tried with the Jacobi method at first but it fails miserably for high density ratios (meaning it converges extremely slowly).

Also, just to be clear, I'm not looking for alternatives (as I have done so quite extensively), I'm trying to debug the issue with the paper in question
Debugging this remotely would deem extremely hard. That's why I was hoping with my initial post to find someone who has actually managed to implement the described algorithm successfully...
XorT is offline   Reply With Quote

Old   June 26, 2020, 10:29
Default
  #52
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
you mean how do i approximate \hat{p} ? as suggested in the paper as \hat{p}^{n+1}=2p^{n}-p^{n-1}
No i mean your approximation for the density gradients.

Quote:
Originally Posted by XorT View Post
I tried with the Jacobi method at first but it fails miserably for high density ratios (meaning it converges extremely slowly).
With two iterations it works quite well for constant density (just tested). But i think you are correct. When the density changes form 1 to 1000, the iterative Poisson solver has to do some work.
__________________
Check out my side project:

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

Old   June 26, 2020, 10:36
Default
  #53
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
No i mean your approximation for the density gradients.

The density gradient is not really needed, but the density at the position where the gradient is definied (between cells). This is done by linear interpolation (simple averaging) of the two neighboring cells as the paper suggests.
XorT is offline   Reply With Quote

Old   June 26, 2020, 11:04
Default
  #54
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
The density gradient is not really needed, but the density at the position where the gradient is definied (between cells). This is done by linear interpolation (simple averaging) of the two neighboring cells as the paper suggests.
Ok i see. From our short conversation i am fairly sure you know what you are doing. Moreover i think you can trust the JCP paper. I think it's more likely to be a bug in your code. I can only recommend that you debug your code very carefully. Perhaps this will help you:

https://en.wikipedia.org/wiki/Rubber_duck_debugging

In the meantime, maybe I or another guy in the forum will come up with another idea.
__________________
Check out my side project:

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

Old   June 26, 2020, 11:05
Default
  #55
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
The density gradient is not really needed, but the density at the position where the gradient is definied (between cells). This is done by linear interpolation (simple averaging) of the two neighboring cells as the paper suggests.



But what about your stability criterion? It seems your problem becomes stiff at high density ratio, I immagine you need some time step constraint to manage that. Furthermore, what about refining the grid?
FMDenaro is offline   Reply With Quote

Old   June 26, 2020, 11:27
Default
  #56
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
Have you also considered this discussion about the interface?
Attached Images
File Type: jpeg 4BEF6DE8-A680-4EBB-89F2-FBB29F43778E.jpeg (32.1 KB, 7 views)
FMDenaro is offline   Reply With Quote

Old   June 30, 2020, 15:08
Default
  #57
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by Eifoehn4 View Post
In the meantime, maybe I or another guy in the forum will come up with another idea.

Yes, maybe someone comes along who did implement this technique successfully. Thank you for your input nevertheless
XorT is offline   Reply With Quote

Old   June 30, 2020, 15:11
Default
  #58
New Member
 
Join Date: Jun 2020
Posts: 27
Rep Power: 5
XorT is on a distinguished road
Quote:
Originally Posted by FMDenaro View Post
But what about your stability criterion? It seems your problem becomes stiff at high density ratio, I immagine you need some time step constraint to manage that. Furthermore, what about refining the grid?

It fluid gains momentum at the density transitions and making the timestep really small seems to just delay the problem. So it doesn't really serve as a solution, unfortunately. I suspect some tiny error in my code, yet i'm unable to find it
XorT is offline   Reply With Quote

Old   June 30, 2020, 15:14
Default
  #59
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
It fluid gains momentum at the density transitions and making the timestep really small seems to just delay the problem. So it doesn't really serve as a solution, unfortunately. I suspect some tiny error in my code, yet i'm unable to find it



Yes, a bug in your code is possible, however reducing the time-step while using a first order upwind should be sufficient to get a stable solution. If that does not happen I would find about a bug for sure.
FMDenaro 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 12:18.