[petsc-users] problem in calling KSPSetComputeRHS and KSPSolve multiple times to solve different RHS.

[Bishesh Khanal]

On Fri, Sep 13, 2013 at 2:47 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

>
> On Sep 13, 2013, at 6:21 AM, Bishesh Khanal <bisheshkh at gmail.com> wrote:
>
> > Dear all,
> > I'm trying to reuse an operator for the linear system multiple times
> with different RHS.
> > But it seems that the RHS is not changed in the KSP object. Here is the
> code to illustrate what I'm doing:
> >
> > Vec mX[3], mB[3];   //Solving Ax=b for 3 different RHS.
> > ...
> > ierr = KSPSetComputeOperators(mKsp,computeMatrix,this);CHKERRQ(ierr);
> > ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);
>
>    Move the line up above the previous line
>
> > ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);
> >
> > for(int i = 0; i<3; ++i) {
> >         mBComponent = i;   //sets different component for RHS.
> >         ierr = KSPSetComputeRHS(mKsp,computeRhs,this);CHKERRQ(ierr);
>
>    Put the above line in the setup code above, not in the loop
>
> >         ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);
>
>    Remove the line above
>
> >         ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
> >         ierr = KSPGetSolution(mKsp,&mX[i]);CHKERRQ(ierr);
> >         ierr = KSPGetRhs(mKsp,&mB[i]);CHKERRQ(ierr);
> >  }
>
>    After you make this changes it should work. If it does not then run
> with -start_in_debugger noxterm and put a break point in computeRHS, make
> sure it is being called.
>
Thanks, I did as you suggested but there is still a problem.
I had kept KSPComputeRHS within the loop so that it gets called to update
the RHS based on the value of mBComponent, using computeRhs function.
Now I kept it outside the loop in the setup code as you suggested and while
doing debugging, I see that the program does go inside the computeRhs
function three times even in this case.
However when I see the results with the following code, it still has the
same rhs values for all three cases.
ierr = PetscObjectSetName((PetscObject)mB[0],"rhsX");CHKERRQ(ierr);
    ierr = VecView(mB[0],viewer);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)mB[1],"rhsY");CHKERRQ(ierr);
    ierr = VecView(mB[1],viewer);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)mB[2],"rhsZ");CHKERRQ(ierr);
    ierr = VecView(mB[2],viewer);CHKERRQ(ierr);

Now I am getting confused with two things:
1) It seems that computeRhs function will be called every time we call
KspSolve. If this is the case, is it the same with KSPSetComputeOperators
too ? Because I do not want the computeMatrix function to be called
repeatedly to assemble the same matrix operator three times. Or there must
be sth I'm missing here.

2) The computeRhs is called three times with different values of
mBComponent, but I can't figure out why the rhs vectors in mB[] array are
same. Here is my computeRhs function:

PetscErrorCode HarmCurvRegSolver::computeRhs(KSP ksp, Vec b, void *ctx)
{
    PetscFunctionBeginUser;
    HarmCurvRegSolver *context = (HarmCurvRegSolver*)ctx;

    PetscErrorCode ierr;
    DM da;
    ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
    DMDALocalInfo info;
    ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);

    PetscScalar ***rhs;
    ierr = DMDAVecGetArray(da,b,&rhs);CHKERRQ(ierr);

    //The RHS should be provided by the model!
    for (PetscInt k = info.zs; k<info.zs+info.zm; ++k) {
        for (PetscInt j = info.ys; j < info.ys+info.ym; ++j) {
            for (PetscInt i = info.xs; i < info.xs+info.xm; ++i) {
                rhs[k][j][i] =
context->getProblem()->getRhsAtPosition(i,j,k,context->getCurrentRhsComponent());

//I've checked well the getRhsAtPosition(...) function! It does give
different values of the 4th argument.
//and getCurrentRhsComponent() just returns the value of mBComponent.
            }
        }
    }

    ierr = DMDAVecRestoreArray(da,b,&rhs);CHKERRQ(ierr);
    ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
    ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

    PetscFunctionReturn(0);

}




>
> >
> > What I expect is to have three different solutions in mX array for the
> different Rhs computed by computeRhs function. However, the result when I
> see using VecView is that the Rhs mB[0], mB[1] and mB[2] are all identical
> (which should not have been the case). And mX also has three identical
> vectors as solution.
> > What am I doing wrong here ?
> >
> >
> > Thanks,
> > Bishesh
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130913/04533dfb/attachment.html>