[FMDenaro]

Originally Posted by XorT

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.

[XorT]

Originally Posted by FMDenaro

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

[FMDenaro]

Originally Posted by XorT

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

[XorT]

Originally Posted by FMDenaro

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.

[FMDenaro]

Originally Posted by XorT

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?

[XorT]

Originally Posted by FMDenaro

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

Originally Posted by FMDenaro

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

[FMDenaro]

Originally Posted by XorT

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

[XorT]

Originally Posted by FMDenaro

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)

[FMDenaro]

Originally Posted by XorT

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.

[Eifoehn4]

Originally Posted by XorT

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:



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

Edit:

Originally Posted by FMDenaro

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

[XorT]

Originally Posted by Eifoehn4

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:



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

you mean how do i approximate ? as suggested in the paper as


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

[Eifoehn4]

Originally Posted by XorT

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.

Originally Posted by XorT

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.

[XorT]

Originally Posted by Eifoehn4

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.

[Eifoehn4]

Originally Posted by XorT

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.

[FMDenaro]

Originally Posted by XorT

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]

Have you also considered this discussion about the interface?

[XorT]

Originally Posted by Eifoehn4

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]

Originally Posted by FMDenaro

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

[FMDenaro]

Originally Posted by XorT

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.