[petsc-users] Split vector with VecGetSubVector

[Matthew Knepley]

On Sun, 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?
>

Run in the debugger with -start_in_debugger, and when it stops type 'where'
and
send all output.

  Thanks,

    Matt


> 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
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141102/f87618c5/attachment.html>