[NLeb]

After solving the momentum predictor MU = -grad(p), the next step in the SIMPLE algorithm is to split MU = AU - H, where A is purely diagonal. This A matrix is then inverted to form the "conductivity" term in the pressure Poisson equation div(A^-1 grad(p)) = div(A^-1 H). To me, this seems like it puts unfair weight on the coordinate axes in which the problem is discretised - the eigenvectors of the A matrix will always correspond with flow along the coordinate axes, even if the flow is actually at 45 deg to the axes. This method does get the correct answer eventually though (with some iteration).

So I wonder: what would happen if we used a block-diagonal A? Could this lead to a reduction in the number of inner correctors / SIMPLE iterations required?

E.g. in 3D, each block would be a full 3x3 matrix, capturing the dependence of the x,y,z-momentum equations on each of the x,y,z velocities at that cell (but not on the neighbours' velocities). The A matrix would still be simple to invert: a full 3x3 block is a bit more expensive to invert than just the diagonal, but this cost is completely insignificant compared to the matrix solve for the pressure correction.

My one concern was that then the A matrix would not be symmetric, leading to an asymmetric pressure matrix, but that could be fixed by simply discarding the skew-symmetric component of A (this is implicitly done if we only allow a diagonal A as usual).

Are there any other issues that could crop up with this idea? Or perhaps there isn't any real benefit to doing this?

I've not yet been able to find any papers or textbooks discussing this - if there are such resources, please could someone point me to them?