Great
> On Jan 26, 2018, at 9:48 AM, Ben Yee <[email protected]> wrote: > > 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. <[email protected]> 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 <[email protected]> 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. <[email protected]> > > 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 <[email protected]> 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
