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

[Jed Brown]

Barry Smith <bsmith at mcs.anl.gov> writes:

>> 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 didn't know it "didn't work" because what we were doing did work.  (I
think we used PCMG interfaces directly; we were mixing unassembled and
assembled representations, plus Galerkin on coarser levels.)

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

Should they be "levels" (numbered?) or named approximations?

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

Presumably.  Is a DM a function space, an operator, or both?

> 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?). 

DMSNES and DMTS effectively provide this, it's just intimately coupled
to the DM.  I advocated before for separating these because function
spaces have meaning independent of operators defined on them, but that's
more of a performance consideration that assumes DM implementations
don't have cheap sharing of common data.

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

SNESSetFunction calls DMSNESSetFunction.

> Or does one every need "lower order operation evaluation"? Is lower
> order linearization enough? 

No, dual-order discretizations are possibly more important for nonlinear
problems.

> 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

One common application of lower order discretizations is in in the
function provided to SNESSetNGS.

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

Yup.

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

I would say quite gradual and with at least a full release cycle of
backward compatibility.  And although I think it is ultimately better
for extensibility, I'm not 100% convinced it is worth the cost to users
to overhaul.  One nice thing about the current interface is that it
forces people to acknowledge that those two arguments could be
different, so they may read or find examples.  If we switch so our
operators have a basket of approximations, lots of users may never learn
about that basket.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 832 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170225/23662de7/attachment.sig>