<div dir="ltr"><div dir="ltr">On Sun, Feb 27, 2022 at 2:36 AM Bojan Niceno <<a href="mailto:bojan.niceno.scientist@gmail.com">bojan.niceno.scientist@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Dear all,<div><br></div><div>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.</div><div><br></div><div>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 <font face="monospace">src/ksp/ksp/tutorials/ex29.c</font>:</div><div><font face="monospace">    MatNullSpace nullspace;<br>    ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);</font></div><div><font face="monospace">    ierr = MatSetNullSpace(J,nullspace);CHKERRQ(ierr);<br>    ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);</font></div><div><br></div><div>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.</div><div><br></div><div>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?</div></div></blockquote><div><br></div><div>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</div><div><br></div><div>   N^T x = 0</div><div><br></div><div>We use this to project out these components at each step of the iterative method. At the end, we get a solution</div><div><br></div><div>  A x = b</div><div><br></div><div>which _also_ satisfies the uniqueness condition</div><div><br></div><div>  N^T x = 0</div><div><br></div><div>When you perturb the matrix, you the the solution to a _different_ linear system</div><div><br></div><div>  (A + sigma I) x = b</div><div><br></div><div>It is close, but not the same, and there is no guarantee (unless you know something about A) that this is close</div><div>to the other solution in a normwise sense.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>    Cheers,</div><div><br></div><div>    Bojan Niceno</div></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>