[petsc-dev] Proposed changes to TS API

[Jed Brown]

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

>> On May 11, 2018, at 6: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.
>> 
>>> 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.
>
>    So this is a backdoor way of allowing the user to provide separate functions for computing the mass matrix and the Jacobian but the current TSSetIJacobian() function only takes one set of matrices in which means that if it returns just the mass matrix it returns it with the storage pattern of the Jacobian, so if the mass matrix is diagonal you get back a sparse matrix with many explicit nonzeros and just the entries on the diagonal. 

Okay, that's messy to work around and relevant to methods like compact
FD where the mass matrix is significantly sparser than the RHS operator.

>    How about the following based on suggestions from several people.
>
>  Keep the current public API TSSetIJacobian() as is
>
>   Add two new public functions:    TSSetMass(Mat, functionmass),  TSSetIJacobianWithOutMass(Mat, functionwithoutmass) (bad name).

What if, when Mass is set, the existing IJacobian is called with shift=0
by any consumers that desire the parts separately?  Then we wouldn't
have to name TSSetIJacobianWithoutMass.  There could be a "don't ever
call me with nonzero shift" option if the IJacobian implementation was
incapable of assembling that part and didn't want to call

  TSGetMassMatrix(ts,&mass);        // honors a reuse policy
  MatAXPY(J,shift,mass,user_knows_best_what_flag_to_use);

>   Provide developer functions    TSComputeMass() which uses functionmass when possible else ijacobian, TSComputeIJacobianWithoutMass() that computes using functionwithoutmass if possible otherwise ijacobian. And change TSComputeIJacobian() so it uses i jacobian if it exists otherwise uses functionmass and functionwithoutmass and combines them automatically. 
>
>   Drawbacks include more complicated internal code (though with the helper functions it will be cleaner especially in the adjoints), power users can provide combined operation, non power users have a simpler interface.