<div dir="ltr"><div dir="ltr">On Sat, Apr 2, 2022 at 8:59 PM Bhargav Subramanya <<a href="mailto:bhargav.subramanya@kaust.edu.sa">bhargav.subramanya@kaust.edu.sa</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 am trying to solve Ax = b in parallel using MUMPS, where x is composed of velocity, pressure, and temperature. There is a null space due to the homogeneous Neumann pressure boundary conditions. </div><div><br></div><div>1. Without explicitly removing the null space, the solution is bounded and robust, and I am able to iterate in time, although I need to check if the solution is physically correct. I am not sure if MUMPS is able to detect the null space in this case. </div></div></blockquote><div><br></div><div>MUMPS should fail with a 0 pivot. I do not understand this.</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>2. However, when I explicitly use MatSetNullSpace as shown below, the solution blows up in 2-3 temporal iterations, resulting in nan. Could you please share your thoughts regarding this problem?</div></div></blockquote><div><br></div><div>Did you check that this is a nullspace vector using <a href="https://petsc.org/main/docs/manualpages/Mat/MatNullSpaceTest.html">https://petsc.org/main/docs/manualpages/Mat/MatNullSpaceTest.html</a> ?</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>Thanks,</div><div>Bhargav</div><div><br></div><div>---------------</div><div>Output for -pc_type svd -pc_svd_monitor </div><div><div>SVD: condition number 3.159257564649e+16, 1 of 1080 singular values are (nearly) zero<br> SVD: smallest singular values: 3.693890048955e-17 2.902795363876e-10 2.902795486255e-10 2.902796682780e-10 7.614221123356e-10<br> SVD: largest singular values : 1.166981958972e+00 1.166981958972e+00 1.166995008014e+00 1.166995008014e+00 1.166995008014e+00</div></div><div>-----------------</div><div>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="color:rgb(0,80,50)">PetscErrorCode</span> <b>SetNullVector</b>(<span style="color:rgb(0,80,50)">Mat</span> MatGlobalLinear,<span style="color:rgb(0,80,50)">Mesh</span> mesh,<span style="color:rgb(0,80,50)">Equations</span> eq){</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="color:rgb(0,80,50)">Vec</span><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>NullVector;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="color:rgb(0,80,50)">PetscInt</span> <span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>numNullVectors,disp,i,j;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="color:rgb(0,80,50)">PetscScalar</span> <span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>constantValue = 1.0,nrm;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="color:rgb(0,80,50)">MatNullSpace</span><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>NullSp = NULL;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo;color:rgb(0,80,50)"><span style="color:rgb(0,0,0)"><span style="white-space:pre-wrap"> </span></span>PetscErrorCode<span style="color:rgb(0,0,0)"><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>ierr;</span></p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo;min-height:19px"><br></p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>PetscFunctionBeginUser;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>numNullVectors = 1;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>disp = mesh.<span style="color:rgb(0,0,192)">numLayers</span>*mesh.<span style="color:rgb(0,0,192)">numEdges</span> + (eq.<span style="color:rgb(0,0,192)">pressure</span>-1)*mesh.<span style="color:rgb(0,0,192)">numCells</span>;</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>MatCreateVecs</b></span>(MatGlobalLinear,&NullVector,NULL);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="color:rgb(127,0,85)"><b>for</b></span>(i=0; i<mesh.<span style="color:rgb(0,0,192)">numLayers</span>; i++){</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span><span style="color:rgb(127,0,85)"><b>for</b></span>(j=0; j<mesh.<span style="color:rgb(0,0,192)">numCells</span>; j++){</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>ierr = VecSetValue(NullVector,j + i*mesh.<span style="color:rgb(0,0,192)">numCells</span> + disp,constantValue,<span style="color:rgb(0,0,192)"><i>INSERT_VALUES</i></span>);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span><span style="white-space:pre-wrap"> </span>}</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>}</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>VecAssemblyBegin</b></span>(NullVector);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>VecAssemblyEnd</b></span>(NullVector);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>VecNorm</b></span>(NullVector,<span style="color:rgb(0,0,192)"><i>NORM_2</i></span>,&nrm);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>VecScale</b></span>(NullVector,1.0/nrm);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo;color:rgb(63,127,95)"><span style="color:rgb(0,0,0)"><span style="white-space:pre-wrap"> </span></span>/* set null space */</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>MatNullSpaceCreate</b></span>(MPI_COMM_WORLD,<span style="color:rgb(0,0,192)"><i>PETSC_FALSE</i></span>,numNullVectors,&NullVector,&NullSp);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>MatSetNullSpace</b></span>(MatGlobalLinear,NullSp);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo;color:rgb(63,127,95)"><span style="color:rgb(0,0,0)"><span style="white-space:pre-wrap"> </span></span>/* clean up */</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>MatNullSpaceDestroy</b></span>(&NullSp);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>ierr = <span style="color:rgb(100,40,128)"><b>VecDestroy</b></span>(&NullVector);CHKERRQ(ierr);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo"><span style="white-space:pre-wrap"> </span>PetscFunctionReturn(0);</p>
<p style="margin:0px;font-variant-numeric:normal;font-variant-east-asian:normal;font-stretch:normal;line-height:normal;font-family:Menlo">}</p></div><div>-----------------<br></div><div><br></div></div>
<br>
<div><hr></div><font face="Arial" size="1">This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.</font></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>