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

[Barry Smith]

On Sep 13, 2013, at 8:35 AM, Bishesh Khanal <bisheshkh at gmail.com> wrote:

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

  Yes

> If this is the case, is it the same with KSPSetComputeOperators too ?

   No, you need to call KSPSetComputeOperators() to trigger a new matrix generation on the next solve.

> 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()); 

put break point here if needed to check that values being set are different each call. Use a tiny da to help debugging.
  

> //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);


    Call VecView() here on b; does it have different values each time or the same?  If it does then you need to run in the debugger and break at the point indicated above and look at the values as they are copied over to make sure they are different each call.

>     ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
>     ierr = VecAssemblyEnd(b);CHKERRQ(ierr);

   You don't need the VecAssemblyBegin/End here.
> 
>     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
> >
> 
>