Efficient schemes for solving the extended Saddle point problem

Question

I am interested in knowing some efficient techniques for solving the following extended Saddle point problem. \begin{align} \begin{bmatrix} A & B^T & C^T \\ B & 0 & 0 \\ C & 0 & 0 \\ \end{bmatrix} \begin{Bmatrix} x \\ y \\ z \\ \end{Bmatrix} = \begin{Bmatrix} F_x \\ F_y \\ F_z \\ \end{Bmatrix} \end{align}

x is velocity, y is pressure, and z is Lagrange multipliers for imposing interface conditions on immersed solids. Such a matrix system is obtained using the fictitious domain/distributed Lagrange multiplier method for incompressible flow problems.

I am looking at two different schemes. In both the schemes, matrices A and B remain the same at each time step but for matrix C the values, as well as the locations of non-zero entries, change from one time step to the other. The number of rows of matrix C is significantly small compared to those of A and B.

Scheme 1: Matrix A is diagonal (and invertible). Since A is diagonal, I am employing a Schur complement based solver using SuperLU solver from the Eigen library. I am wondering if the efficiency can be improved further since only matrix C, which is significantly smaller in size compared to B, changes at each time step.

Scheme 2: Matrix A is sparse and symmetric (and invertible). At the moment I use the sparse direct solver PARDISO. This is certainly not an efficient approach since matrices A and B remain the same. I do not know what type of Schur complement approach would be efficient.

I am looking for some alternative strategies that can be efficient. I very much appreciate any inputs in this regard.

Note: I tried Uzawa type schemes in the past. They are neither efficient nor robust for generic problems; I do not like to fiddle with the relaxation parameters in the Uzawa scheme and its variants.

Answer 1

For both schemes: I am going to lump $B$ and $C$ into a single matrix $B$ for convenience of communicating the idea. Likewise, $y$ and $z$ are now $y$, and $F_y$ and $F_z$ are $F_y$.

$Ax+B^Ty=F_x$ or $x = A^{-1}(F_x-B^Ty)$

$Bx=F_y$ or $BA^{-1}(F_x-B^Ty)=F_y$

$BA^{-1}B^Ty=BA^{-1}F_x-F_y$

Note that $BA^{-1}B^T$ is symmetric and invertible, and can be factorized using a cholesky decomposition. To construct $BA^{-1}B^T$ (assuming the $A$ is amenable to a cholesky decomposition):

$$BA^{-1}B^T = BL^{-T}L^{-1}B^T=(L^{-1}B^T)^T(L^{-1}B^T)$$ which is amenable to a symmetric Rank-K update, which is fast, efficient, and present in most linear algebra libraries. In BLAS, it is usually called something like syrk; in Eigen, it is called by invoking selfAdjointView<Lower>().rankUpdate(). If $A$ is diagonal, this is much simpler.

EDIT: Regarding your comment that $B$ is actually quite large:

$$ \begin{pmatrix} B \\ C \\ \end{pmatrix}A^{-1} \begin{pmatrix} B^T C^T \\ \end{pmatrix}= \begin{pmatrix} BA^{-1}B^T & BA^{-1}C^T \\ CA^{-1}B^T & CA^{-1}C^T \\ \end{pmatrix} $$

The top left block of this matrix will be constant. The top right is redundant (in a lower triangular storage scheme). If you use syrk to update the top left and bottom right corners, they will store and compute only the lower triangular portion, which saves a lot of time and space.

You can either accept the resulting block matrix and solve it, or you can use a block inversion scheme to determine the action of that inverse using the cholesky factorizations of the top left and bottom right corners. This may be useful since you only have to factorize the top left corner once.

Answer 2

Disclaimer: I wrote this answer in a rush, may not be up to the standards. Also, I don't have enough information about your problem to give more detailed advice. However, the information presented here should be enough to get you going.

I will suggest you to use MINRES with an appropriate preconditioner.

For a moment, let's assume that this is not a twofold saddle point problem, and you have the following problem to solve

$\begin{equation} \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} U \\ P \end{bmatrix} = \begin{bmatrix} F \\ G \end{bmatrix}. \end{equation}$

I am assuming that you are using an inf-sup stable discretization of the Navier-Stokes equations, (for scheme 2) which would imply that $BA^{-1}B^T$ and $Q_p$, the mass matrix defined on the pressure space, are spectrally equivalent. Moreover, you can prove that $Q_p$ and $diag(Q_p)$ are spectrally equivalent, and spectral equivalence is a transitive relation. Also, with scheme 2, $A$ is just a vector Laplacian, hence you could use (algebraic) multigrid to very efficiently invert it. Such a preconditioner can be described through the action of its inverse as

$P^{-1} = \begin{bmatrix} AMG(A) & 0 \\ B & diag(Q_p)^{-1} \end{bmatrix}.$

Note that I am using $P^{-1}$ to mean as the action of the inverse of

$P= \begin{bmatrix} A & 0 \\ B & diag(Q_p) \end{bmatrix}.$

If $diag(Q_p)$ is too slow in terms of iterations, you can use $AMG(Q_p)$ to obtain a better preconditioner (in terms of iterations) at slightly higher computational cost.

Now, we need to deal with the second saddle point problem. I am not familiar with immersed solid literature, but I think there must be an inf-sup condition for that part too. Then by using Theorem 3.1 of (Howell JS, Walkington NJ (2011) Inf-sup conditions for twofold saddle point problems.), you can find another spectrally equivalent matrix -potentially a mass matrix depending on problem parameters-, let's call it $Z$. Finally, the following preconditioner probably would be optimal for your problem

$P= \begin{bmatrix} A & 0 & 0\\ B & diag(Q_p) & 0\\ C & 0 & Z \end{bmatrix}.$

Depending on the form of $Z$, you may be able to approximate it with its diagonal part or using an algebraic multigrid method (I suggest BoomerAMG). Or maybe if the nonzero structure of $Z$ is nice, for example block-diagonal with small size diagonal blocks, you can use a direct solver at low computational cost.