[petsc-users] Split vector with VecGetSubVector

[Barry Smith]

  You may also want to consider using the PETSc PCFIELDSPLIT preconditioner; this is designed for solving problems like Stokes and many others. With PCFIELDSPLIT you build one matrix that defines the entire problem over the velocities and pressures together and then it internally builds up solvers for the blocks so you do not need to monkey with all those details yourself.

  Barry

> On Nov 2, 2014, at 7:19 AM, Massimiliano Leoni <leoni.massimiliano1 at gmail.com> wrote:
> 
> 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