[WIP] [linear] Add first draft of SIMPLE preconditioner

SIMPLE

Goal: solve

\begin{pmatrix}
A & B\\
C & D
\end{pmatrix}
\begin{pmatrix}
u\\
p
\end{pmatrix}
=
\begin{pmatrix}
f\\
g
\end{pmatrix}

Equations

Exact velocity: ~~~~~~u = u^* + \delta u ~~~~~~~ (1)

Exact pressure: ~~~~~~p = p^* + \delta p ~~~~~~~ (2)

Exact velocity equation: ~~~~~~A u + B p = f ~~~~~~~ (3)

Approximate velocity equation using some pressure guess p^*: ~~~~~~A u^* + B p^* = f ~~~~~~~ (4)

Substract: (3-4) : ~~~~~~A (u - u^*) + B (p-p^*) = 0 ~~~~~~~ (5)

~~~~~~A \delta u + B \delta p = 0 ~~~~~~~ (6)

Assumption (SIMPLE): ~~~~~~diag(A) \approx A ~~~~~~~ (7)

Rearrange (6) and insert (7): ~~~~~~\delta u = - diag(A)^{-1} B \delta p~~~~~~~ (7)

Mass conservation: ~~~~~~C u = g~~~~~~~ (8)

~~~~~~C (u^* + \delta u) = g~~~~~~~ (9)

~~~~~~C (u^* + - C diag(A)^{-1} B \delta p) = g~~~~~~~ (10)

~~~~~~C u^* = g + C diag(A)^{-1} B \delta p~~~~~~~ (11)

Schur complement ~~~~~~C diag(A)^{-1} B \delta p = C u^* - g ~~~~~~~ (12)

Algorithm

1.) Solve (4) for u^*

2.) Solve (12) for \delta p

3.) Get \delta u from (7)

4.) update p=p^* + \alpha \delta p

5.) update u= u^* +\delta u

Naive implementation of https://www.cs.umd.edu/~elman/papers/tax.pdf

Converges, but is slower than Uzawa.

The internal GMRES solver to invert the Schur complement is not preconditioned yet (Richardson, w = 1). Richardson seems to be the only preconditioner that does not require an assembled matrix.

Maybe we can write our own Jacobi preconditioner? If one can apply the action of the Schur complement without a matrix (what we do right now), is it also possible to apply the inverse of its diagonal without having the matrix?

Udpate : There seem to be various matrix-free preconditioning methods

http://www2.cs.cas.cz/~tuma/ps/cutu06.pdf

https://hal.inria.fr/inria-00074080/document

http://www2.cs.cas.cz/semincm/lectures/2009-05-15-DuintjerTebbens.pdf

https://arxiv.org/pdf/1805.11930.pdf

https://arxiv.org/pdf/2006.06052.pdf

https://gitlab.dune-project.org/exadune/exadune-appl/-/blob/feature/blockiterative-preconditioner-FlexibleGMRES/src/blockiterativesolver/blockdiagonalpreconditioner.hh

A very simple (and probably expensive) option to test if preconditioning helps inverting the approximate Schur \tilde{S} complement would be

\tilde{S}^{-1}_{0,0}  = \left(\tilde{S} \cdot (1, 0, ... ,n)^T \right) \cdot (1, 0, ... ,n)^T \\

\tilde{S}^{-1}_{1,1}  = \left(\tilde{S} \cdot (0, 1, ... ,n)^T \right) \cdot (0, 1, ... ,n)^T \\
....

With that, one could construct a Jacobi preconditioner after doing n matrix vector operations.