<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Feb 5, 2015 at 4: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:0 0 0 .8ex;border-left:1px #ccc 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></blockquote><div><br></div><div>What the person who told you this may have meant is that CG/GMRES will</div><div>not fail with your singular operator. However, you can get any one of the </div><div>infinite possible solutions depending on your rhs vector. Set the nullspace.</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
[*] 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="HOEnZb"><font color="#888888">Massimiliano<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">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></div>