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

[Ben Yee]

Hi Barry,

This seems to work for me.  I defined the following function in our code
repository, and I call it every time I use KSPSetOperators to update the
operator for my KSP object with PCMG:

PetscErrorCode PCMGSoftReset(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; 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);
    }
  }
  pc->setupcalled = PETSC_FALSE;
  PetscFunctionReturn(0);
}

(Note that I had to add "pc->setupcalled = PETSC_FALSE".)  Thanks for the
help!

On Wed, Jan 24, 2018 at 12:19 PM, Smith, Barry F. <bsmith at mcs.anl.gov>
wrote:

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


-- 
Ben Yee

NERS PhD Candidate, University of Michigan
B.S. Mech. & Nuclear Eng., U.C. Berkeley
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180126/8f76e561/attachment-0001.html>