[petsc-users] Changing the operator for PCMG while leaving interpolation matrices unchanged

[Smith, Barry F.]

   This is off the top of my head and my not be useful. Copy PCReset_MG() which is 

PetscErrorCode PCReset_MG(PC pc)
{
  PC_MG          *mg        = (PC_MG*)pc->data;
  PC_MG_Levels   **mglevels = mg->levels;
  PetscErrorCode ierr;
  PetscInt       i,n;

  PetscFunctionBegin;
  if (mglevels) {
    n = mglevels[0]->levels;
    for (i=0; i<n-1; i++) {
      ierr = VecDestroy(&mglevels[i+1]->r);CHKERRQ(ierr);
      ierr = VecDestroy(&mglevels[i]->b);CHKERRQ(ierr);
      ierr = VecDestroy(&mglevels[i]->x);CHKERRQ(ierr);
      ierr = MatDestroy(&mglevels[i+1]->restrct);CHKERRQ(ierr);
      ierr = MatDestroy(&mglevels[i+1]->interpolate);CHKERRQ(ierr);
      ierr = VecDestroy(&mglevels[i+1]->rscale);CHKERRQ(ierr);
    }

    for (i=0; i<n; i++) {
      ierr = MatDestroy(&mglevels[i]->A);CHKERRQ(ierr);
      if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
        ierr = KSPReset(mglevels[i]->smoothd);CHKERRQ(ierr);
      }
      ierr = KSPReset(mglevels[i]->smoothu);CHKERRQ(ierr);
    }
  }
  PetscFunctionReturn(0);
}

And change it to only destroy mats on each level ierr = MatDestroy(&mglevels[i]->A);CHKERRQ(ierr); Keep everything else.

Not that you will need to include #include <petsc/private/pcmgimpl.h>  in your code.

> On Jan 24, 2018, at 9:53 AM, Ben Yee <bcyee at umich.edu> wrote:
> 
> Yes, the nonzero structure changes in a way that cannot easily be predicted.  I could overallocate, add some extra entries to the first matrix, and reuse that matrix, but I would have to roughly double the number of allocated entries, leading to an undesirable increase in memory.  I think that one alternative would be to just reset the KSP object when I change the operator, but, if I'm not mistaken, this would destroy all the interpolation operators and I'd have to redefine all of them.
> 
> On Wed, Jan 24, 2018 at 10:41 AM, Smith, Barry F. <bsmith at mcs.anl.gov> wrote:
> 
>   Ben
> 
>     This is a good question, I'll take a look at how difficult it would be to add this feature.
> 
>      I take it the nonzero structure of the outer matrix changes over time?
> 
> 
>    Barry
> 
> 
> > On Jan 24, 2018, at 9:12 AM, Ben Yee <bcyee at umich.edu> wrote:
> >
> > Hi,
> >
> > I have multiple matrices (with different nonzero structures) for which I would like to use the same KSP operator (KSPRICHARDSON with PCMG, Galerkin process for generating coarse grid operator).  I want to use the same PCMG structure for all the matrices (same interpolation, same restriction, same number of levels, same smoothing techniques, etc.), so I tried using KSPSetOperators to just change the operator.  My interpolation matrices are MPIAIJ matrices that I have defined myself and provided to PCMG.  My hope was that, with the new operator I provided, PETSc would recompute all of the coarse/intermediate grid operators using the Galerkin process with the old interpolation matrices and the new fine-grid operator.
> >
> > However, I get the wrong solution after changing the operator with KSPSetOperators.  After some extensive digging around the PETSc source code, it seems that the operators in mg.c are not updated by KSPSetOperators.  When I do a MatView on mglevels->A or the operator in mglevels->smoothu or mglevels->smoothd in PCMGMCycle_Private() of mg.c, it has the matrix values of the old matrix.  From what I understand, the only place where mglevels->A (used to compute the residual) is updated is in PCApplyRIchardson_MG() of mg.c, from the following lines of code:
> >  74:   for
> >  (i=0; i<levels; i++) {
> >
> >  75:     if
> >  (!mglevels[i]->A) {
> >
> >  76:       KSPGetOperators
> > (mglevels[i]->smoothu,&mglevels[i]->A,NULL);
> >
> >  77:       PetscObjectReference((PetscObject
> > )mglevels[i]->A);
> >
> >  78:
> >     }
> >
> >  79:   }
> > It seems to me that the pointer mglevels[i]->A is only set once, when it is undefined at the beginning of the iteration scheme.  I think there are more issues as well, but this is the easiest one to point out.  (For example, it doesn't seem like it is possible to update the operator for mglevels[i]->smoothd from KSPSetOperators on the outermost KSP object, and, if I try to set it manually via PCMGGetSmoother, I get errors if the nonzero structures of the coarse-grid operators have changed, which tends to be the case when the fine-grid operator's nonzero structure changes).
> >
> > Is there something I'm doing wrong?  Is there a way to do this so that I can reuse PCMG with a new operator?
> >
> > Thanks!
> >
> > --
> > Ben Yee
> >
> > NERS PhD Candidate, University of Michigan
> > B.S. Mech. & Nuclear Eng., U.C. Berkeley
> 
> 
> 
> 
> -- 
> Ben Yee
> 
> NERS PhD Candidate, University of Michigan
> B.S. Mech. & Nuclear Eng., U.C. Berkeley