Solving saddle point problem having non-invertible top-left block with a PETSc nested matrix
My system is a symmetric FE problem with lagrange multipliers:
$Z=\begin{pmatrix}A & C^T \ C & 0\end{pmatrix}$
The matrix $A$ is positive semi-definite, non-invertible. The whole matrix is invertible.
I am working on the development of a finite element software where I want to solve such problems with PETSc in parallel. At this time I would like the $Z$ matrix to be a MATNEST in which the $A$ and $C$ matrices are assembled independently.
With the representative simple code at the end of this message (1D laplacian with Lagrange multipliers), and using a PCFIELDSPLIT preconditioner, the solver diverges. If I artificially force the matrix $A$ to be invertible (adding the command-line option -force_invertible), then I get the right solution.
So, what can I do to solve such a system with a nested matrix?
Another question: is it possible to use direct methods (Mumps) with nested matrices in parallel?
#include "petscsys.h" /* framework routines */
#include "petscvec.h" /* vectors */
#include "petscmat.h" /* matrices */
#include "petscksp.h"
#include <vector>
#include <string>
#include <iostream>
#include <numeric>
// Try to solve with Petsc a simple 1D laplacian problem (or a series of springs)
// o-////-o-////-o-////-o-////-o-////-o-////-o-////-o-////-o ...
//
// The boundary conditions (Dirichlet) are imposed using Lagrange multipliers.
// The system to be solved is of the form:
//
// Z x = y
// --^-- -^- -^-
// [A Ct] [u] = [b]
// [C 0 ] [l] [d]
//
// A is positive semi-definite (1 eigenvalue is 0)
// Z is indefinite invertible
static char help[] = "Saddle point problem: scalar 1D laplacian with lagrange multipliers.\n\n";
int main(int argc,char **argv)
{
// Initialization ------------------------------------------------ //
MPI_Init(NULL, NULL);
PetscErrorCode ierr;
ierr = PetscInitialize(&argc,&argv,NULL,help);CHKERRQ(ierr);
int rank, nproc;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &nproc);
PetscViewer viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE);
PetscBool forceInvertible=PETSC_FALSE;
PetscOptionsGetBool(NULL,NULL, "-force_invertible", &forceInvertible, NULL);
// Problem definition and dof spliting among processes ----------- //
double k=1.;
std<IP_ADDRESS>vector<int> globalDofs={0, 1, 2, 3, 4, 5, 6, 7, 8};
std<IP_ADDRESS>vector<std<IP_ADDRESS>pair<int, double>> globalDirichlet={{0, 0.},
{8, 1.}}; // first: global dof number, second: imposed value
auto start=[&](int rank){
int q=(globalDofs.size()-1)/nproc;
int r=(globalDofs.size()-1)%nproc;
return q*rank + std<IP_ADDRESS>min(rank, r);
};
std<IP_ADDRESS>vector<int> dofs;
for (int i=start(rank); i<=start(rank+1); ++i) {
dofs.push_back(globalDofs[i]);
}
std<IP_ADDRESS>vector<int> isLocal(dofs.size(), 1);
if (rank! 935