Preconditioning improves the conditioning of the Krylov operator.
\[ \begin{gather*} (P^{-1} A) x = P^{-1} b \\ \{ P^{-1} b, (P^{-1}A) P^{-1} b, (P^{-1}A)^2 P^{-1} b, \dotsc \} \end{gather*} \]
\[ \begin{gather*} (A P^{-1}) P x = b \\ \{ b, (P^{-1}A)b, (P^{-1}A)^2b, \dotsc \} \end{gather*} \]
Various preconditioning strategies are presented in
Feel++ abstracts the PETSc library and provides a subset (sufficient in most cases) to the PETSc features. It interfaces with the following PETSc libraries: Mat
, Vec
, KSP
, PC
, SNES
.
Vec:
Vector handling libraryMat:
Matrix handling libraryKSP:
Krylov SubSpace library implements various iterative solversPC:
Preconditioner library implements various preconditioning strategiesSNES:
Nonlinear solver library implements various nonlinear solve strategiesAll linear algebra are encapsulated within backend which interfaces PETSc (petsc
), Eigen sparse (eigen
) or dense (eigen_dense
)
The Default backend is petsc
.
We start with the quickstart Laplacian example, recall that we wish to, given a domain \(\Omega\), find \(u\) such that
\[ -\nabla \cdot (k \nabla u) = f \mbox{ in } \Omega \subset \mathbb{R}^{2}, u = g \mbox{ on } \partial \Omega \]
shell> mpirun -np 2 feelpp_qs_laplacian --ksp-monitor=1 --ksp-view=1 0 KSP Residual norm 8.953261456448e-01 1 KSP Residual norm 7.204431786960e-16 KSP Object: 2 MPI processes type: gmres GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement GMRES: happy breakdown tolerance 1e-30 maximum iterations=1000 tolerances: relative=1e-13, absolute=1e-50, divergence=100000 left preconditioning using nonzero initial guess using PRECONDITIONED norm type for convergence test PC Object: 2 MPI processes type: shell Shell: linear system matrix = precond matrix: Matrix Object: 2 MPI processes type: mpiaij rows=525, cols=525 total: nonzeros=5727, allocated nonzeros=5727 total number of mallocs used during MatSetValues calls =0 not using I-node (on process 0) routines
We now turn to the quickstart Stokes example, recall that we wish to, given a domain \(\Omega\), find \((\mathbf{u},p) \) such that
\[ -\Delta \mathbf{u} + \nabla p = \mathbf{ f},\, \nabla \cdot \mathbf{u} = 0 \mbox{ in } \Omega, \mathbf{u} = \mathbf{g} \mbox{ on } \partial \Omega \]
This problem is indefinite, Matrices.
\[ \label{eq:2} M = \begin{pmatrix} A & B\\ B^T & 0 \end{pmatrix} = \begin{pmatrix} I & 0\\ B^T C & I \end{pmatrix} \begin{pmatrix} A & 0\\ 0 & - B^T A^{-1} B \end{pmatrix} \begin{pmatrix} I & A^{-1} B\\ 0 & I \end{pmatrix} \]
\[ \label{eq:3} P_1 = \begin{pmatrix} \tilde{A}^{-1} & B\\ B^T & 0 \end{pmatrix}, \quad P_2 = \begin{pmatrix} \tilde{A}^{-1} & 0\\ 0 & \tilde{S} \end{pmatrix},\quad P_3 = \begin{pmatrix} \tilde{A}^{-1} & B\\ 0 & \tilde{S} \end{pmatrix} \]
where \(\tilde{S} \approx S^{-1} = B^T A^{-1} B\) and \(\tilde{A}^{-1} \approx A^{-1}\)