[petsc-dev] Proposed changes to TS API

[Matthew Knepley]

On Fri, May 11, 2018 at 7:05 PM, Jed Brown <jed at jedbrown.org> wrote:

> "Smith, Barry F." <bsmith at mcs.anl.gov> writes:
>
> >    Here is MY summary of the discussion so far.
> >
> > 1) the IFunction/IJacobian interface has its supporters. There is valid
> argument that for certain cases it can be more efficient than the proposed
> alternative; but this seems largely a theoretical believe at this time
> since there are no comparisons between nontrivial (or even trivial) codes
> that use the assembly directly the mass shift plus Jacobian vs the assembly
> separately and MatAXPY the two parts together.  As far as I can see this
> potential performance is the ONLY benefit to keeping the
> IFunction/IJacobian() interface? Please list additional benefits if there
> are any?
>
> PCShell
>
> > 2) The IFunction/IJacobian approach makes it impossible for TS to access
> the mass matrix alone.
>
> Without wasted work, yes.
>
> > 3) But one can access the IJacobian matrix alone by passing a shift of
> zero
> >
> > 4) TSComputeIJacobian() is an ugly piece of shit code that has a
> confusing name since it also can incorporates the RHS Jacobian.
>
> If you get rid of that, then every implicit integrator will need to
> handle the RHSFunction/RHSJacobian itself.  It will be significantly
> more code to maintain.
>
> > 4a) the TSComputeIJacobian() has issues for linear problems because it
> overwrites the user provided Jacobian and has hacks to deal with it
>
> Yes, however that choice reduces memory usage which is sometimes an
> important factor.
>
> > 5) There is no standard way to solve M u = F() explicitly from the
> IFunction() form, (and cannot even with expressed with the explicit RHS
> approach) the user must manage the M solve themselves inside their
> RHSFunction.
>
> We could write this as an IMEX method with IFunction providing M u and
> RHSFunction providing F.  I think this could be specialized for constant
> M and provided automatically for any explicit method.


We need this for PyLith, and are doing what you suggest by hand right now.
Automating it would
go a long way toward removing objection to the current division I think.

  Matt


>
> > 6) There is some support for adding two new function callbacks that a)
> compute the mass matrix and b) compute the Jacobian part of the IJacobian.
> This appears particularly useful for implementing adjoints.
> >
> > 7) If we added the two new interfaces the internals of an already
> > overly complicated TS become even more convoluted and unmaintainable.
>
> We could split the shift into ashift and bshift (elided as always 1.0
> now), then users could opt to skip work if one of those was nonzero.
> That would be a modest generalization from what we have now and would
> enable any desired optimization.  Integrators that wish to take
> advantage of M not changing could interpret a flag saying that it was
> constant and then always call the IJacobian with ashift=0.  That would
> be unintrusive for existing users and still enable all optimizations.
> It's only one additional parameter and then optimized user code could
> look like
>
>   for (each element) {
>     Ke[m][m] = {};
>     if (ashift != 0) {
>       Ke += ashift * (mass stuff);
>     }
>     if (bshift != 0) {
>       Ke += bshift * (non-mass stuff);
>     }
>   }
>
> Naive user code would elide the conditionals, but would still work with
> all integrators.
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20180512/d1b7b712/attachment.html>