[petsc-users] Split vector with VecGetSubVector

[Massimiliano Leoni]

@Matt: I luckily managed to solve the problem myself, I started reading 
through the PETSc manual and re-wrote cleaner code. It runs, but yields a 
wrong result.

@Barry: this definitely looks interesting! I am on a university assignment and 
I'll often have to solve this kind of problems, so this PCFIELDSPLIT looks 
just like what is needed. Also, people at FEniCS dev are very interested in 
it, so I want to go deeper.

I read the documentation and the tutorial and I wrote this version, starting 
from my original working code and the tutorial:
void StokesIdBlTrPrecSegregated::solve(PETScVector& x, const PETScVector& b)
{
    // PCFIELDSPLIT version
    // "the whole world is talking about it,
    // and I've got to know why"

	[inits, bcs, and other stuff I'll move to the constructor asap]


//    KSP            *subksp;
//    PetscInt       n = 1;
    IS isu;
    int Nu = Mv->size(0);
    ISCreateStride(PETSC_COMM_WORLD,Nu,0,1,&isu);
    IS isp;
    int Np = Mp->size(0);
    ISCreateStride(PETSC_COMM_WORLD,Np,0,1,&isp);

    PC pc;

    KSP ksp;
    KSPCreate(PETSC_COMM_WORLD,&ksp);
    KSPGetPC(ksp, &pc);
    PCFieldSplitSetIS(pc, "0", isu);
    PCFieldSplitSetIS(pc, "1", isp);

/*    if (s->userPC) {			//not yet clear what's this for
      PCFieldSplitSchurPrecondition(pc, PC_FIELDSPLIT_SCHUR_PRE_USER, s->myS);
    }
    if (s->userKSP) {
      PCSetUp(pc);
      PCFieldSplitGetSubKSP(pc, &n, &subksp);
      KSPSetOperators(subksp[1], s->myS, s->myS, SAME_PRECONDITIONER);
      PetscFree(subksp);
    }
*/
    Mat Pr(P->mat());

    KSPSetOperators(ksp,Pr,Pr,DIFFERENT_NONZERO_PATTERN);
    KSPSetFromOptions(ksp);
    KSPSolve(ksp,b.vec(),x.vec());

    KSPDestroy(&ksp);
}


Now, this runs and yields the correct solution, but how can I be sure it is 
actually doing splitting blocks? This looks rather incomplete to me.
What's missing in this code?
More importantly, where can I read more on this PCFIELDSPLIT?

Thanks for the precious help,
Massimiliano

P.s. Shall I change the thread subject to another? Or start another thread?



In data domenica 2 novembre 2014 15:30:26, hai scritto:
>   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/ex
> > 70.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