[petsc-users] Split vector with VecGetSubVector

[Massimiliano Leoni]

Hi everyone,
  this is my first email so I'll use this sentence to say hello from Milan, 
Italy :)

I am very new to PETSc, I always "used" it through another library, FEniCS, 
but now I need to use its interface directly.

Straight to my problem, I have this working code:

// PETScVector is a FEniCS container for a PETSc Vec, which can be retrieved 
// with vec() method; similarly P->mat().
void solve(PETScVector& x, const PETScVector& b)
{
    Mat Pr(P->mat());

    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD,&ksp);
    KSPSetOperators(ksp,Pr,Pr,DIFFERENT_NONZERO_PATTERN);
    KSPSolve(ksp,b.vec(),x.vec());

    KSPDestroy(&ksp);
}

which simply solves with a linear system.

Now I want to do the same, but separating the system [the matrix is block 
diagonal] into two parts [velocity and pressure, that's Stokes!].

My try, basing on my basic understanding of the documentation

void solve(PETScVector& x, const PETScVector& b)
{
    // I understand I need these Index Sets to split a vector
    IS isu;
    int Nu = Mv->size(0);		// number of components of velocity
    ISCreateStride(PETSC_COMM_SELF,Nu,0,1,&isu);

    Vec buv;		// to contain velocity rhs
    VecCreate(PETSC_COMM_WORLD,&buv);	// I am not sure when to use 
    VecSetSizes(buv,PETSC_DECIDE,Nu);		// VecCreate and VecSetSizes
    VecGetSubVector(b.vec(),isu,&buv);

    Mat Mvp(Mv->mat());			// again, retrieve a Mat
    Vec xuv;		// to contain velocity solution
    VecCreate(PETSC_COMM_WORLD,&xuv);
    VecSetSizes(xuv,PETSC_DECIDE,Nu);

    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD,&ksp);
    KSPSetOperators(ksp,Mvp,Mvp,DIFFERENT_NONZERO_PATTERN);
    KSPSolve(ksp,buv,xuv);		[** segfault! **]

    KSPDestroy(&ksp);

    [pressure part is analogous] 


    Vec finalVec;		// to re-assemble whole vector
    VecCreate(PETSC_COMM_WORLD,&finalVec);
    VecSetSizes(finalVec,PETSC_DECIDE,Nu+Np);

    IS iss[2] = {isu,isp};
    Vec vecs[2] = {xuv,xpv};
    VecCreateNest(PETSC_COMM_SELF, 2, iss, vecs, &finalVec);

    VecCopy(finalVec,x.vec());
}


The problem is that it segfaults at KSPSolve.
Since the structure is not too far from my working code above, I would guess I 
am missing something when extracting the subvectors, maybe some initializing 
procedure?

I read the tutorial code http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex70.c.html in which this is used, but 
it's a little too complicated for me now.

I am using petsc 3.4.2 [debian testing].

Thanks for help!
Massimiliano