[petsc-dev] DMKSP and pcfieldsplit

[Lawrence Mitchell]

Hi Matt, all,

> On 21 Feb 2019, at 18:17, Matthew Knepley <knepley at gmail.com> wrote:
> 
> Here is why I did this. I wanted this to work
> 
>   https://bitbucket.org/petsc/petsc/src/4dbc1805575afffed4e440f1353fcfccbc893081/src/snes/examples/tutorials/ex62.c#lines-1062
> 
> The DMKSP is initialized by SNES to have a function which says "get your matrix from SNES". Without this, the subproblem could
> not make the matrix. Now I changed the description of exactly what matrix it is making and stored that in the subDM.
>  
> What is the right way to think about fixing this? I think that if the DM can create a field decomposition, it should also be able to transfer user callbacks from the parent to the subdms.
> 
> That way, the DM (IOW me) can control how that happens.
> 
> Yes, that is what I thought I was doing above. Preserving the callbacks that had been set in the original DM.
> Maybe clarify what you think should happen.
>  
> Or do I have the way this is supposed to work all wrong. Is the ComputeOperators callback I set meant to be completely agnostic, and it should get all the problem-specific information from the DM?
> 
> I am doing this, but it does not have to be the only way things work. What I do not understand is how this change breaks your stuff.
> If you were setting the DMKSP compute callback, it should just override what I am doing above. If not, why does me setting it here
> screw it up?

I've done some more digging, and the setup that breaks is:

SNES with pcfieldsplit inside and rediscretised multigrid inside the fieldsplit.

What is going on is, I think, the following:

SNESSolve sets up a ComputeOperators on the KSP:

    KSP ksp;
    ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
    ierr = KSPSetComputeOperators(ksp,KSPComputeOperators_SNES,snes);CHKERRQ(ierr);
    ierr = DMCoarsenHookAdd(snes->dm,DMCoarsenHook_SNESVecSol,DMRestrictHook_SNESVecSol,snes);CHKERRQ(ierr);

How does this callback know where to linearise around? Well, KSPComputeOperators_SNES does this:

  if (dmsave == snes->dm) X = snes->vec_sol; /* We are on the finest level */
  else {                                     /* We are on a coarser level, this vec was initialized using a DM restrict hook */
    ierr = DMGetNamedGlobalVector(snes->dm,"SNESVecSol",&Xnamed);CHKERRQ(ierr);


And the DMCoarsenHookAdd and DMRestrictHook_SNESVecSol ensure that /if we coarsen the DM associated with the SNES/ that state vector is transferred to the coarser levels.

However, when we do PCFieldsplit and DMCreateFieldDecomposition, this restriction hook and connection to the state vector is gone.

As a result, when I rediscretise on the coarse grids, I'm always linearising around a zero state.

So I suppose that what should happen is that somewhere in fieldsplit, we need to transfer SNES state over and setup the appropriate restrict hooks.

Cheers,

Lawrence