[arkie87]

Originally Posted by sbaffini

In contrast, a simple Symmetric Gauss Seidel (just plain SGS, not as smoother or preconditioner), works ridiculously well in the coupled density based solver.

I once found myself at a dinner with a group of prominent linear algebra professors, except that I didn't know it (I knew they were somehow involved in math/engineering, but didn't know exactly on which topic), and then I suddenly dropped this bomb (using plain SGS in a production code that solves large system of equations)... they were all ... I was all "I know, but it just works"... it was fun

I assume this is because advection term is dominant, and as long as you sweep downstream it will work extremely well? I assume advection term is dominant because density based solver is usually used in high mach/compressible flows? Or would it work well for natural convection?

[sbaffini]

Originally Posted by arkie87

I assume this is because advection term is dominant, and as long as you sweep downstream it will work extremely well? I assume advection term is dominant because density based solver is usually used in high mach/compressible flows? Or would it work well for natural convection?

I didn't mention it but, of course, GMRES+ILU0 works better, basically always. Where better is less outer iterations for fixed inner iterations (or even less inner iterations).

The unexpected was that SGS basically always works, in a way or another, and the difference with GMRES+ILU0 is not order of magnitudes, but O(1), maybe 10 sometimes. Yet SGS is really faster, so the overall time is down to comparable in most cases.

It is so for both convection and diffusion dominated problems, and the main explanation I found is that, no matter what, you are still solving a sytem linearized at a given pseudo-time/outer iteration. If it is, say, the first outer iteration, there is little to gain from oversolving it. And it is the same for most of the initial phase. See, for example, here:

https://aisberg.unibg.it/retrieve/ha...ntinuation.pdf
http://www.ii.uib.no/~trond/publicat...rs/Inexact.pdf
https://users.wpi.edu/~walker/Papers...1996,16-32.pdf

In practice, you mostly need just 1-2 orders in the residual reductions per outer iteration, and SGS can actually give them to you. SGS will fail in going any further, but that is not required in most cases.

[arkie87]

Originally Posted by sbaffini

No, forget about inverting the whole matrix at full level, which is feasible only in MATLAB for very few cells.

Ok, forgotten.

Originally Posted by sbaffini

Real CFD codes actually use some form of iterative algebraic solver to solve the full system (either for the n segregated or the 1 coupled).

Yeah, my SIMPLE code used GMRES, though most of my problems were small enough that I just used direct solver.

Originally Posted by sbaffini

In such solvers you typically need to divide the residual or the rhs by the diagonal coefficient aii (think about gauss-seidel). Now, what happens in a coupled solver is that your diagonal coefficient is not anymore a single scalar, but the nxn matrix Aii, so instead of dividing the rhs by aii you actually multiply it by the inverse of Aii, that's where the inversion comes in for coupled solvers. But Aii blocks (and Aij as well, for that matter) are already full, and so are their inverse. So, to summarize, you only invert the diagonal nxn blocks, you do it once per outer iteration, and you actually store the inverse in place, so the overall overhead is relatively limited.

Makes sense.

Originally Posted by sbaffini

Imagine a segregated solver, that solves for u, then v, then p. I make an example for finite volumes because it is more natural for me, but the same apply to other methods. In solving for u, you make a loop over all the faces of the grid, retrieve the variables on the left and right of a face and compute the fluxes which are then assigned to the left and right residuals, or something along these lines. Then you do your loop over the cells (as opposed to faces) to compute additional source terms, etc.
Now, no matter how well ordered are the faces and cells, there is always some cache issue in retrieving variables from cells and putting back some other values. That is, at some point (unless the problem is very small), you will need a cell value that is out of cache, during each of those loops. Now, solving a segregated problem implies that those problems are encountered n times. Also, you kind of end up retrieving and computing the same variable multiple times (imagine the diffusion coefficients, for example).

With a coupled solver you do only one face loop and one cell loop, retrieve all the variables only once, compute what you need once, and put back stuff once. All the variables of a cell, in my ordering, are stored together and don't incur cache misses when retrieved together. Actually, this effectively reduces the number of cells in cache (because there are more variables per cell), but the cache misses are reduced.

In theory, a segregated solver could do the same i.e. retrieve the info on all equations for cell ij instead of just for equation k. Or do I misunderstand?

Originally Posted by sbaffini

Yes, it seems so, but the u-v coupling in the paper would be in the term mf, which is in afuu and afvv (see eq. 23). They decide to skip it, but including it would have added no cost. Note that with Rhie-Chow mf also contains the pressure gradient (see equations 24 and 25 in the paper).

Not sure I understand you here. Eq. (23) shows afuv = afvu=0 not afuu=afvv=0

Originally Posted by sbaffini

Yes, I guess that coupling equations in a staggered framework comes with its own problems.

Yeah, i will probably fully read the paper and see if it is worth redoing my code to be collocated. For now, I am not using blocks, and I can solve the system using GMRES if i replace the continuity equation with the pressure equation from SIMPLER algorithm (in order to have a non-zero diagonal).

Originally Posted by sbaffini

As said above, inversion is only involved for the diagonal blocks. Coupling with neighbor cells is the same as for segregated solvers.

I am still trying to understand how this works. you multiply block Aii by the inverse of Aii, I assume you also multiply Aij (across, not down) by the inverse of Aii for all Aij=/=0? It is sort of like combining U V and P into one variable and performing Guass seidel on that?

[arkie87]

Originally Posted by sbaffini

I didn't mention it but, of course, GMRES+ILU0 works better, basically always. Where better is less outer iterations for fixed inner iterations (or even less inner iterations).

The unexpected was that SGS basically always works, in a way or another, and the difference with GMRES+ILU0 is not order of magnitudes, but O(1), maybe 10 sometimes. Yet SGS is really faster, so the overall time is down to comparable in most cases.

It is so for both convection and diffusion dominated problems, and the main explanation I found is that, no matter what, you are still solving a sytem linearized at a given pseudo-time/outer iteration. If it is, say, the first outer iteration, there is little to gain from oversolving it. And it is the same for most of the initial phase. See, for example, here:

https://aisberg.unibg.it/retrieve/ha...ntinuation.pdf
http://www.ii.uib.no/~trond/publicat...rs/Inexact.pdf
https://users.wpi.edu/~walker/Papers...1996,16-32.pdf

In practice, you mostly need just 1-2 orders in the residual reductions per outer iteration, and SGS can actually give them to you. SGS will fail in going any further, but that is not required in most cases.

Ok, I guess i mistook and didnt realize SGS is different than GS. From this reference https://math.stackexchange.com/quest...idel-iteration it seems like SGS inverts the (D+L) matrix to get x_k+1/2, and then uses x_k+1/2 to get x_k by inverting (D+U) matrix.

It is sort of like assuming that the east and north values are fixed and solving for east/west, and then reversing that assumption. Because it uses a matrix inverse, i expect:

(1) That low frequency errors will be dissipated fast (unlike regular GS)
(2) That it is much less expensive to get the inverse of (D+L) and (D+U) than (D+L+U). I suppose (D+L) is like an advection only problem i.e. it is parabolic not elliptic, so a sweep in the right direction will solve it even with regular GS.

[sbaffini]

Originally Posted by arkie87

Not sure I understand you here. Eq. (23) shows afuv = afvu=0 not afuu=afvv=0

afuv and afvu are 0 exactly because the inter-velocity coupling is in mf which only multiplies the current variable. Said otherwise, u-equation sees v only trough mf=mf(u,v), but mf only multiplies u in the u equation (there is no other v except in mf).

Originally Posted by arkie87

It is sort of like combining U V and P into one variable and performing Guass seidel on that?

YES. Not sure for other methods, but for GS it is EXACTLY THAT.

[sbaffini]

Originally Posted by arkie87

Ok, I guess i mistook and didnt realize SGS is different than GS. From this reference https://math.stackexchange.com/quest...idel-iteration it seems like SGS inverts the (D+L) matrix to get x_k+1/2, and then uses x_k+1/2 to get x_k by inverting (D+U) matrix.

It is sort of like assuming that the east and north values are fixed and solving for east/west, and then reversing that assumption. Because it uses a matrix inverse, i expect:

(1) That low frequency errors will be dissipated fast (unlike regular GS)
(2) That it is much less expensive to get the inverse of (D+L) and (D+U) than (D+L+U). I suppose (D+L) is like an advection only problem i.e. it is parabolic not elliptic, so a sweep in the right direction will solve it even with regular GS.

The only difference is that you sweep back and forth instead of always using the same sweep direction

[arkie87]

Originally Posted by sbaffini

afuv and afvu are 0 exactly because the inter-velocity coupling is in mf which only multiplies the current variable. Said otherwise, u-equation sees v only trough mf=mf(u,v), but mf only multiplies u in the u equation (there is no other v except in mf).

Right. It is explicit though, rather than implicit.

[arkie87]

Originally Posted by sbaffini

The only difference is that you sweep back and forth instead of always using the same sweep direction

I guess there is no need to invert the matrix then. Inverting the matrix i assume is slower than simply using a for loop over all the rows of the matrix.

But in that case, i would expect it to be terrible at dampening out the low frequency errors.

[sbaffini]

Originally Posted by arkie87

I guess there is no need to invert the matrix then. Inverting the matrix i assume is slower than simply using a for loop over all the rows of the matrix.

But in that case, i would expect it to be terrible at dampening out the low frequency errors.

That's SGS

[arkie87]

Originally Posted by sbaffini

That's SGS

Thanks again for all your help.

I am looking into Rhie Chow and collocated grid. I have the CFD textbook by Ferziger, and it describes a term like Rhie Chow that is proportional to the third derivative of the pressure. It also appears in a similar form as Rhie Chow but is not labeled as Rhie Chow, so I am not quite sure if it is the same thing or not. For the purpose of this post, I am going to assume it is the same thing, but kindly please let me know if it is not.

As far as i can tell, the basic idea behind the Rhie Chow method is to solve the momentum equations using the same CVs as the continuity and other scalars, but add corrections to momentum and/or continuity equations for consistency and energy conservation and to prevent oscillating pressure fields. That said I have a few questions, if you dont mind:

(1) At this point, I dont know how the cell-centered pressure gradient in the momentum equations is computed-- is it simply computed using central difference [(P_E-P_W)/(2*dx)], which results in coefficients of P_P = 0 (since no P_P term appears in the expression above), or is there a Rhie chow correction term (and if so, how is it defined)?

(2) From Ferziger, it seems like there is a Rhie Chow correction in the momentum equation at the eastern face (which is probably used for continuity), but i dont know if there is any correction at the cell center?

(3) The velocity corrections are applied to the faces, seemingly without a Rhie Chow correction. The only difference is that the A_p coefficients are computed by interpolation onto the faces rather than getting the A_p coefficient using the values of face velocities (as in a staggered grid). I assume then that the cell velocity corrections are applied by averaging the face corrections.

[arjun]

Originally Posted by sbaffini

In contrast, a simple Symmetric Gauss Seidel (just plain SGS, not as smoother or preconditioner), works ridiculously well in the coupled density based solver.

I once found myself at a dinner with a group of prominent linear algebra professors, except that I didn't know it (I knew they were somehow involved in math/engineering, but didn't know exactly on which topic), and then I suddenly dropped this bomb (using plain SGS in a production code that solves large system of equations)... they were all ... I was all "I know, but it just works"... it was fun

This is indeed true. I actually first came to know about it when I was starting off my career and we had to use a code written at Tohoku university. They use this for high speed simulation.