[XorT]

Originally Posted by Eifoehn4

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?

Hmm, I'm affraid I don't know enough about it to answer this. However, what I did verify is that after applying the pressure correction, the divergence of the velocity field is in fact zero everywhere. Regardless of this, though, the velocities grow indefinetely every timestep.

[Eifoehn4]

Originally Posted by FMDenaro

To be honest, the paper is not very clear...

That's true.

[XorT]

Originally Posted by FMDenaro

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?

I think you are mistaken. The constant poisson equation is solved directly and accurately. The solution exactly solves the equation. Also the velocity field has zero divergence after the correction.


PS.: yes, i assume h=1 for simplicity

[FMDenaro]

Originally Posted by XorT

I think you are mistaken. The constant poisson equation is solved directly and accurately. The solution exactly solves the equation. Also the velocity field has zero divergence after the correction.


PS.: yes, i assume h=1 for simplicity

I just considered the equations in the paper. Try yourself by hand, build the discrete div v on the staggered grid and check the result. You will see that the discrete 5-point FD Poisson equation appears. But if you use the spectral solution that is not verified exactly.

[XorT]

Originally Posted by FMDenaro

I just considered the equations in the paper. Try yourself by hand, build the discrete div v on the staggered grid and check the result. You will see that the discrete 5-point FD Poisson equation appears. But if you use the spectral solution that is not verified exactly.

I don't know what you mean exactly. Are you saying that the use of the poisson equation is wrong ? Because the kind of solver used to solve it should be completely irrelevant. Either it solves the equation or it doesn't (in which case it is not really a solver, but merely an approximation of some kind).
I'm not getting why we're stuck on the FFT based solver topic.

[FMDenaro]

Originally Posted by XorT

I don't know what you mean exactly. Are you saying that the use of the poisson equation is wrong ? Because the kind of solver used to solve it should be completely irrelevant. Either it solves the equation or it doesn't (in which case it is not really a solver, but merely an approximation of some kind).
I'm not getting why we're stuck on the FFT based solver topic.

Maybe you are not aware of the basic principle about approximate and exact projection methods. One of the fundamental is the use of the same discrete operators in each step. Using the FFT means you are adopting a spectral interpolation for the pressure. Then the pressure gradient should be computed from the spectral polynomial not from the discret gradient shown in Eqs. (39) and (40). That is not clear to me, how can be ensured the divergence-free field up to machine precision?

Again, use (39) and (40) to build the Poisson equation. You get a discrete FD equation while the paper addressed the spectral interpolation for the continuous Poisson problem.

[XorT]

Originally Posted by FMDenaro

Maybe you are not aware of the basic principle about approximate and exact projection methods. One of the fundamental is the use of the same discrete operators in each step. Using the FFT means you are adopting a spectral interpolation for the pressure. Then the pressure gradient should be computed from the spectral polynomial not from the discret gradient shown in Eqs. (39) and (40). That is not clear to me, how can be ensured the divergence-free field up to machine precision?

Again, use (39) and (40) to build the Poisson equation. You get a discrete FD equation while the paper addressed the spectral interpolation for the continuous Poisson problem.

Isn't the whole paper using the discretized version ? Isn't the discretized poisson equation solved exactly (by the solver of your choice). I don't see this dependent on a spectral solver at all. Where do you see this statement ? You can use whatever solver you want as long as it solves the fully discretized poisson equation. They even employ MG methods as an alternative.
I really don't get it...

[FMDenaro]

Originally Posted by XorT

Isn't the whole paper using the discretized version ? Isn't the discretized poisson equation solved exactly (by the solver of your choice). I don't see this dependent on a spectral solver at all. Where do you see this statement ? You can use whatever solver you want as long as it solves the fully discretized poisson equation. They even employ MG methods as an alternative.
I really don't get it...

The FFT is not just a solver! It implies a functional representation of the function p by means of trigonometric functions. If you start from Lap p that is a continuous representation of the operator, equivalent to Div Grad. In projection methods you have to discretize Grad and Div in a congruent way but that is not necessarily equivalent to the discrete Lap operator.

That is the basic theory of the projection methods you can find in the original Chorin's paper as well as in the textbooks of CFD.
The paper addressed (39) and (40) as if the Div Grad operator was discretized by second order central formula. But, actually, it seem that Lap p is directly solved by the FFT.

[XorT]

Originally Posted by FMDenaro

The FFT is not just a solver! It implies a functional representation of the function p by means of trigonometric functions. If you start from Lap p that is a continuous representation of the operator, equivalent to Div Grad. In projection methods you have to discretize Grad and Div in a congruent way but that is not necessarily equivalent to the discrete Lap operator.

That is the basic theory of the projection methods you can find in the original Chorin's paper as well as in the textbooks of CFD.
The paper addressed (39) and (40) as if the Div Grad operator was discretized by second order central formula. But, actually, it seem that Lap p is directly solved by the FFT.

Well, I'm not expert on this, but as I said, I verified the results. Here is some proof for two test-cases (once I hard transition, once a random field for the RHS). Computed the pressure, computed the gradient, computed the 5-point laplacian and computed the sum of squared differences to the RHS as error (see above laplacian image). The error is in the real of machine precision I'd say, at least it is so tiny that this is not the problem here, I'm quite sure. (These are the DST versions for Dirichlet BC, btw)



https://ibb.co/LtQHKcY
https://ibb.co/HVwdCBF

[FMDenaro]

Originally Posted by XorT

Well, I'm not expert on this, but as I said, I verified the results. Here is some proof for two test-cases (once I hard transition, once a random field for the RHS). Computed the pressure, computed the gradient, computed the 5-point laplacian and computed the sum of squared differences to the RHS as error (see above laplacian image). The error is in the real of machine precision I'd say, at least it is so tiny that this is not the problem here, I'm quite sure. (These are the DST versions for Dirichlet BC, btw)



https://ibb.co/LtQHKcY
https://ibb.co/HVwdCBF

From your figures I don't see anything I am talking about. The problem is not the solution of the Laplacian but the computation of div v^n+1 in each node.

Could you post the field of the result ?
Again, this is a mathematical conclusion from Eq.(39) and (40). Take the discrete divergence of the LHS and set to zero. Then take the same divergence of the RHS and develop the discrete pressure equation. This is the FD equation that should be solved, not the spectral representation of the laplacian of the pressure. The article is not clear and I do not understand what is done.

[FMDenaro]

Maybe the key is in the source term, it is modified to include the FD laplacian contribution?

[XorT]

Well, we seem to run in circles.
Maybe, for reference, check the standard approach to pressure correction used in fluid dynamics for computer graphics here: http://users.encs.concordia.ca/~grog...cs/fluid-5.pdf (section 5.3.3). They do exactly the same thing. Compute the divergence (cell centered) from the staggered velocities (at cell boundaries), use the divergence (at cell center) to build the RHS of the poisson equation which is then solved to yield pressure (at cell center). Then the gradient of the pressure (which is at cell boundary) is applied to the staggered velocities (also at cell boundary).

Are you saying this a problem by itself already ?

[FMDenaro]

I am stating that you cannot mix different discrete operatore in these steps!

[Eifoehn4]

Originally Posted by FMDenaro

I am stating that you cannot mix different discrete operatore in these steps!

Maybe we should look at this more from the point of view of a computer-graphic scientist.

Originally Posted by XorT

Well, I'm not expert on this, but as I said, I verified the results. Here is some proof for two test-cases (once I hard transition, once a random field for the RHS). Computed the pressure, computed the gradient, computed the 5-point laplacian and computed the sum of squared differences to the RHS as error (see above laplacian image). The error is in the real of machine precision I'd say, at least it is so tiny that this is not the problem here, I'm quite sure. (These are the DST versions for Dirichlet BC, btw)

https://ibb.co/LtQHKcY
https://ibb.co/HVwdCBF

Did you also perform a convergence test, e.g. with a manufactured solutions? This would be a first step to ensure that you at least converge to the right solution no matter what order. Moreover what did you do so far to tackle the instability issue?

[FMDenaro]

Originally Posted by Eifoehn4

Maybe we should look at this more from the point of view of a computer-graphic scientist.



Did you also perform a convergence test, e.g. with a manufactured solutions? This would be a first step to ensure that you at least converge to the right solution no matter what order.

But that jcp paper is about a CFD topic... I would assess if the finally velocity is discretely divergence-free at machine accuracy or at local truncation error accuracy.
From a point of a computer graphics I can accept also a first order stable solution but as here the problem is the instability ...

[Eifoehn4]

Originally Posted by FMDenaro

But that jcp paper is about a CFD topic... I would assess if the finally velocity is discretely divergence-free at machine accuracy or at local truncation error accuracy.
From a point of a computer graphics I can accept also a first order stable solution but as here the problem is the instability ...

Yes that's true. Such inconsistencies may also influence the stability.

[XorT]

Originally Posted by FMDenaro

I am stating that you cannot mix different discrete operatore in these steps!

so http://users.encs.concordia.ca/~grog...cs/fluid-5.pdf (section 5.3.3) does it wrong as well ? because this is well established standard way of doing it when it comes to CFD in CG.

Originally Posted by FMDenaro

But that jcp paper is about a CFD topic... I would assess if the finally velocity is discretely divergence-free at machine accuracy or at local truncation error accuracy.

Both papers are about a CFD topic. And as I said before, I did assess that the final velocity is discretely divergence-free at machine accuracy.

[XorT]

Originally Posted by Eifoehn4

Did you also perform a convergence test, e.g. with a manufactured solutions? This would be a first step to ensure that you at least converge to the right solution no matter what order. Moreover what did you do so far to tackle the instability issue?

I'm trying to simplify the simulation and scenario to home in on the source of the error numerically. I apply no forces, use really small timesteps, but this does not prevent the blow-up UNLESS i use constant density (in which case everything works as expected). The latter also being the reason to doubt, that this discussion is going in the right direction. It is clear to me that the problem is in the "hack" of approximating the pressure for the variable density from past pressures and density gradients. Not what we're discussing about here for days now

[FMDenaro]

Originally Posted by XorT

so http://users.encs.concordia.ca/~grog...cs/fluid-5.pdf (section 5.3.3) does it wrong as well ? because this is well established standard way of doing it when it comes to CFD in CG.






Both papers are about a CFD topic. And as I said before, I did assess that the final velocity is discretely divergence-free at machine accuracy.

Look, exactly what I am stating, the pressure equation is not solved by a spectra interpolation

[XorT]

Originally Posted by FMDenaro

Look, exactly what I am stating, the pressure equation is not solved by a spectra interpolation

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