[petsc-dev] [petsc-users] Multigrid with defect correction

[Barry Smith]

> On Feb 25, 2017, at 7:35 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> 
> Barry Smith <bsmith at mcs.anl.gov> writes:
>>     Do you think this is a reasonable approach or am I missing
>>     something fundamental? I am assuming generally for the "higher
>>     order" DM the Mat it returns is a MATSHELL or a new matrix class
>>     built on "tensor contractions" and that kind of nonsense. I don't
>>     want to do all the coding and then have it turn out that it is
>>     totally useless for CEED etc.
> 
> Well, this is exactly what we do in pTatin.

  So why didn't you fix it in PETSc at that time?

> 
> I would ask, why just two discretizations?  I've always thought a better
> interface would be for Mats (and any other operators/functions) to have
> optional approximations or supplementary data attached.  We can do this
> with PetscObjectCompose, but that's hard to work with in a structured
> way.  Anyway, I think I would rather just have the Amat with ability to
> attach one or more Pmats.

  Understood. I don't like the API as "attach one or more Pmats to a Mat" but instead prefer to think of it as a linear operator that has one or more "levels of approximation" that can optionally be indicated; which one you wish to use is up to the algorithm using the linear operator. Changing this in PETSc would require some conceptional changes and refactoring in PETSc. Would simplify things like SNESSetJacobian() because we would not need to track mat and pmat. But would not be terribly difficult for the PETSc developers or user.

  If we used this model, instead of adding a new DM to the interface do we keep the one DM but have a way to "attach" other "lower order" DMs to the that one DM that can be accessed? 

DMCreateMatrix(dm,&mat);
pmat = mat;
DMGetLowerDM(dm,&ldm); 
if (ldm) {
  DMCreateMatrix(ldm,&pmat);
}
KSP/SNES/TSSetIJacobian(ts,mat,pmat,...); 

more generally DMGetLowerDM(dm,PetscOrderIndicater i,&ldm);  Similarly MatGetLowerMat(mat,PetscOrderIndicater i,&lmat)?

----- an aside
  BTW: this introduces a possible "missing" object in PETSc. We have as a fundamental object the Mat which is a linear operator, we do not have a fundamental object for a "general" operator, instead we accept user provided functions directly (SNESSetFunction(), TSSetIFunction() etc).  I've resisted the PetscOperator because I felt it was overly abstract (and thus confusing to most users) and is more difficult to implement cleanly since the operator may depend on more than one input vector (TSSetIFunction depends on u and udot; would we have a PetscOperator2 that had two input vectors?). I did try the PF object but it is very limited and rarely used. Certainly I think for a wide class of users, especially Fortran, SNESSetFunction() is understandable but 

PetscOperatorCreate(&operator)
PetscOperatorSetRawFunction(operator, function())
SNESSetOperator(snes,operator)

becomes gibberish and so we would lose many potential users. Of course we could introduce the PetscOperator and maintain the simplicity of the current interface by having things like SNESSetFunction() automatically handle the conversion to the "operator" interface and have a new SNESSetOperator() for the general case. I hate having this kind of "wrapper" code because it usually bit decays over time and I don't like have multiple equivalent interfaces for the same effect (how do you document them, how does the user decide which one to use? I hate perl :-), there should be one way to do each thing).

Or does one every need "lower order operation evaluation"? Is lower order linearization enough? Currently with mat and pmat we only support lower order linearization, we don't for example allow passing a higher order and lower order nonlinear function both to SNESSetFunction. If the lower order operator is never needed, only the lower order linearization then PetscOperator (which would allow managing several orders of general operator) would not be needed so API for users is simpler. Or we could introduce lower order linear operators now and not worry about lower order function evaluations for years?
--- end of the aside

I also don't disagree with Matt, should we just "complete" support for two levels of approximation by having proper two DMs for TS/KSP/SNES, see how that goes rather than jumping directly from a (far from perfect since users need to manually do ugly stuff since we only have one DM) two level support (as we currently have) to multilevel support? I won't change the current interface to do this instead I would just add TS/SNES/KSPSetDMs()

Imagine poor Ed's book if we suddenly change stuff like this in PETSc :-(  So the question is what should/could we do now to improve flexibility without requiring Ed to start all over again? That is, how gradually we need to introduce this kind of overhaul, and with how much backward compatibility?


Barry