<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Sun, Nov 2, 2014 at 7:19 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 everyone,<br>
  this is my first email so I'll use this sentence to say hello from Milan,<br>
Italy :)<br>
<br>
I am very new to PETSc, I always "used" it through another library, FEniCS,<br>
but now I need to use its interface directly.<br>
<br>
Straight to my problem, I have this working code:<br>
<br>
// PETScVector is a FEniCS container for a PETSc Vec, which can be retrieved<br>
// with vec() method; similarly P->mat().<br>
void solve(PETScVector& x, const PETScVector& b)<br>
{<br>
    Mat Pr(P->mat());<br>
<br>
    KSP ksp;<br>
    KSPCreate(PETSC_COMM_WORLD,&ksp);<br>
    KSPSetOperators(ksp,Pr,Pr,DIFFERENT_NONZERO_PATTERN);<br>
    KSPSolve(ksp,b.vec(),x.vec());<br>
<br>
    KSPDestroy(&ksp);<br>
}<br>
<br>
which simply solves with a linear system.<br>
<br>
Now I want to do the same, but separating the system [the matrix is block<br>
diagonal] into two parts [velocity and pressure, that's Stokes!].<br>
<br>
My try, basing on my basic understanding of the documentation<br>
<br>
void solve(PETScVector& x, const PETScVector& b)<br>
{<br>
    // I understand I need these Index Sets to split a vector<br>
    IS isu;<br>
    int Nu = Mv->size(0);               // number of components of velocity<br>
    ISCreateStride(PETSC_COMM_SELF,Nu,0,1,&isu);<br>
<br>
    Vec buv;            // to contain velocity rhs<br>
    VecCreate(PETSC_COMM_WORLD,&buv);   // I am not sure when to use<br>
    VecSetSizes(buv,PETSC_DECIDE,Nu);           // VecCreate and VecSetSizes<br>
    VecGetSubVector(b.vec(),isu,&buv);<br>
<br>
    Mat Mvp(Mv->mat());                 // again, retrieve a Mat<br>
    Vec xuv;            // to contain velocity solution<br>
    VecCreate(PETSC_COMM_WORLD,&xuv);<br>
    VecSetSizes(xuv,PETSC_DECIDE,Nu);<br>
<br>
    KSP ksp;<br>
    KSPCreate(PETSC_COMM_WORLD,&ksp);<br>
    KSPSetOperators(ksp,Mvp,Mvp,DIFFERENT_NONZERO_PATTERN);<br>
    KSPSolve(ksp,buv,xuv);              [** segfault! **]<br>
<br>
    KSPDestroy(&ksp);<br>
<br>
    [pressure part is analogous]<br>
<br>
<br>
    Vec finalVec;               // to re-assemble whole vector<br>
    VecCreate(PETSC_COMM_WORLD,&finalVec);<br>
    VecSetSizes(finalVec,PETSC_DECIDE,Nu+Np);<br>
<br>
    IS iss[2] = {isu,isp};<br>
    Vec vecs[2] = {xuv,xpv};<br>
    VecCreateNest(PETSC_COMM_SELF, 2, iss, vecs, &finalVec);<br>
<br>
    VecCopy(finalVec,x.vec());<br>
}<br>
<br>
<br>
The problem is that it segfaults at KSPSolve.<br>
Since the structure is not too far from my working code above, I would guess I<br>
am missing something when extracting the subvectors, maybe some initializing<br>
procedure?<br></blockquote><div><br></div><div>Run in the debugger with -start_in_debugger, and when it stops type 'where' and</div><div>send all output.</div><div><br></div><div>  Thanks,</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">
I read the tutorial code <a href="http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex70.c.html" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex70.c.html</a> in which this is used, but<br>
it's a little too complicated for me now.<br>
<br>
I am using petsc 3.4.2 [debian testing].<br>
<br>
Thanks for help!<br>
<span class="HOEnZb"><font color="#888888">Massimiliano<br>
</font></span></blockquote></div><br><br clear="all"><div><br></div>-- <br>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>