In both cases, it is like you are solving a nonsingular system with a matrix B. With MatNullSpace, B=A-e*e' where e=ones(n,1) normalized, and with your approach it is B=A+sigma*I with sigma=1e-6. The first approach shifts the zero eigenvalue, while in the second approach all eigenvalues are shifted. Jose

El 27 feb 2022, a las 8:36, Bojan Niceno <[email protected]> escribió:

Dear all,

I have coupled PETSc with my computational fluid dynamics (CFD) solver for incompressible flows where the most computationally intensive part is a solution of the linear system for pressure - which is singular.

A simple call to PETSc solvers resulted in divergence, as expected, but things work when I set the null space for the pressure matrix as demonstrated in src/ksp/ksp/tutorials/ex29.c: MatNullSpace nullspace; ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr); ierr = MatSetNullSpace(J,nullspace);CHKERRQ(ierr); ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);

However, the effect of setting the null space as described above, has almost the same effect (convergence history is almost the same) as if when I multiply each diagonal of the system matrix with (1.0 + 1.0e-6), i.e., desingularize the matrix by making it slightly diagonally dominant.

I prefer the former solution as the latter one seems a bit like an ad-hoc patch and I am not sure how general it is, but I wonder, from a mathematical point of view, is it the same thing? Any thoughts on that?

Cheers,

Bojan Niceno

On Sun, Feb 27, 2022 at 2:36 AM Bojan Niceno < [email protected]> wrote:

Dear all,

I have coupled PETSc with my computational fluid dynamics (CFD) solver for incompressible flows where the most computationally intensive part is a solution of the linear system for pressure - which is singular.

A simple call to PETSc solvers resulted in divergence, as expected, but things work when I set the null space for the pressure matrix as demonstrated in src/ksp/ksp/tutorials/ex29.c: MatNullSpace nullspace; ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr); ierr = MatSetNullSpace(J,nullspace);CHKERRQ(ierr); ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);

However, the effect of setting the null space as described above, has almost the same effect (convergence history is almost the same) as if when I multiply each diagonal of the system matrix with (1.0 + 1.0e-6), i.e., desingularize the matrix by making it slightly diagonally dominant.

I prefer the former solution as the latter one seems a bit like an ad-hoc patch and I am not sure how general it is, but I wonder, from a mathematical point of view, is it the same thing? Any thoughts on that?

I will give a slightly different explanation than Jose. When you set a nullspace N, it tells us what space the solution must be orthogonal to N^T x = 0 We use this to project out these components at each step of the iterative method. At the end, we get a solution A x = b which _also_ satisfies the uniqueness condition N^T x = 0 When you perturb the matrix, you the the solution to a _different_ linear system (A + sigma I) x = b It is close, but not the same, and there is no guarantee (unless you know something about A) that this is close to the other solution in a normwise sense. Thanks, Matt

Cheers,

Bojan Niceno

-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>