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