<div dir="ltr"><div class="gmail_default" style="font-family:arial,helvetica,sans-serif">If your linear system comes from a Laplace equation with Neumann boundary conditions (singular matrix) it means for a given b there are infinite possible x.<br></div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif"><br></div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif">So, you get one x and calculate the corresponding b. When you solve Ax=b, you might not get the exact same x, since there are infinite possible x for the same b.</div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif"><br></div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif"><br></div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif"><br></div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif"><br></div><div class="gmail_default" style="font-family:arial,helvetica,sans-serif"><br></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Feb 5, 2015 at 8:14 AM, Massimiliano Leoni <span dir="ltr"><<a href="mailto:leoni.massimiliano1@gmail.com" target="_blank">leoni.massimiliano1@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Hi petsc-users,<br>
I stumbled across a curious issue while trying to setup some solver for<br>
Navier-Stokes.<br>
In one sentence, I tested a solver computing b = Ax with x known, then solving<br>
the system with b to find the original x back, but it's not working.<br>
<br>
In detail, I want to implement a matrix free method for a matrix C2 = D^T<br>
M_L^{-1} D, where M_L^{-1} is the inverse lumped of M; here's what I do<br>
<br>
[I'm sorry for posting code, but I suspect the error could be quite subtle...]<br>
<br>
////// define utility struct<br>
typedef struct _C2ctx<br>
{<br>
Mat D;<br>
Vec Mdiag;<br>
} C2ctx;<br>
<br>
////// define multiplication<br>
int C2mult(Mat A,Vec x,Vec y)<br>
{<br>
C2ctx* ctx;<br>
MatShellGetContext(A,&ctx);<br>
Mat D = ctx->D;<br>
Vec Mdiag = ctx->Mdiag;<br>
int Nu;<br>
VecGetSize(Mdiag,&Nu);<br>
<br>
Vec tmp;<br>
VecCreate(PETSC_COMM_WORLD,&tmp);<br>
VecSetType(tmp,VECSTANDARD);<br>
VecSetSizes(tmp,PETSC_DECIDE,Nu);<br>
<br>
MatMultTranspose(D,x,tmp);<br>
VecPointwiseDivide(tmp,tmp,Mdiag);<br>
MatMult(D,tmp,y);<br>
<br>
VecDestroy(&tmp);<br>
return 0;<br>
}<br>
<br>
Later on, I assemble D and M, and then<br>
<br>
////// lump M into _Mdiag<br>
Vec _Mdiag;<br>
VecCreate(PETSC_COMM_WORLD,&_Mdiag);<br>
VecSetType(_Mdiag,VECSTANDARD);<br>
VecSetSizes(_Mdiag,PETSC_DECIDE,Nu);<br>
MatGetRowSum(_M,_Mdiag);<br>
<br>
////// set context<br>
C2ctx ctx;<br>
ctx.D = _D;<br>
ctx.Mdiag = _Mdiag;<br>
<br>
////// create _C2 [= A]<br>
Mat _C2;<br>
MatCreateShell(PETSC_COMM_WORLD,Np,Np,PETSC_DETERMINE,PETSC_DETERMINE,&ctx,&_C2);<br>
MatShellSetOperation(_C2,MATOP_MULT,(void(*)(void))C2mult);<br>
<br>
////// create _b2 [= b]<br>
Vec _b2;<br>
VecCreate(PETSC_COMM_WORLD,&_b2);<br>
VecSetType(_b2,VECSTANDARD);<br>
VecSetSizes(_b2,PETSC_DECIDE,Np);<br>
<br>
////// create _deltaP [= x]<br>
Vec _deltaP;<br>
VecCreate(PETSC_COMM_WORLD,&_deltaP);<br>
VecSetType(_deltaP,VECSTANDARD);<br>
VecSetSizes(_deltaP,PETSC_DECIDE,Np);<br>
<br>
////// set _deltaP = 1 and _b2 = _C2*_deltaP, then change _deltaP for check<br>
VecSet(_deltaP,1.0);<br>
MatMult(_C2,_deltaP,_b2);<br>
VecSet(_deltaP,2.0);<br>
<br>
////// setup KSP<br>
KSP _ksp2;<br>
KSPCreate(PETSC_COMM_WORLD,&_ksp2);<br>
KSPSetOperators(_ksp2,_C2,_C2,DIFFERENT_NONZERO_PATTERN);<br>
KSPSetType(_ksp2,"cg");<br>
<br>
KSPSolve(_ksp2,_b2,_deltaP);<br>
<br>
According to my understanding, now _deltaP should be a nice row on 1s, but<br>
what happens is that it has some [apparently] random zeros here and there.<br>
<br>
Additional infos:<br>
[*] the problem comes from the cavity problem with navier-stokes. This<br>
particular step is to solve for the pressure and the matrix _C2 is singular<br>
[it's a laplacian] with nullspace the constants. I was told I can ignore this<br>
and avoid setting the nullspace as long as I use a CG or GMRES solver.<br>
<br>
[*] the number of "random zero entries" is very high with P2-P1 elements, and<br>
is significantly lower with P1+Bubble - P1.<br>
<br>
[*] Matrices M and D are automatically assembled by FEniCS, not by me.<br>
<br>
Can anybody please advice on what I am doing wrong? What could be causing<br>
this?<br>
<br>
Thanks in advance<br>
<span class=""><font color="#888888">Massimiliano<br>
</font></span></blockquote></div><br></div></div>