"Smith, Barry F." <bsm...@mcs.anl.gov> writes:
>> On May 11, 2018, at 6:05 PM, Jed Brown <j...@jedbrown.org> wrote:
>>
>> "Smith, Barry F." <bsm...@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
>
> Please explain what this means, I have no clue.
If you are using matrix-free preconditioning, SNESGS, and anything else
that really looks into the problem structure, then the
IFunction/IJacobian interface is much cleaner than having separate
functions for mass and non-mass parts and needing the solver to learn
about it by interrogating a different data structure like MatShell.
>>> 2) The IFunction/IJacobian approach makes it impossible for TS to access
>>> the mass matrix alone.
>>
>> Without wasted work, yes.
>
> The computation of two Jacobians correct?
Yeah, though it could be cached if that was a performance bottleneck.
>>> 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.
>
> I don't propose "getting rid of it", I suggest perhaps the code could be
> refactored (I don't know how) to have something more streamlined.
Okay.
>>> 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.
>
> Then also provide new TS developer routines such as
> TSComputeMass(), TSComputeJacobianWithoutMass() (bad name) so that
> TS code would look cleaner and not have a bunch of ugly calls to
> TSComputeIJacobian with shifts of 1 and 0 around that I suspect
> Hong doesn't like.
Either way; I don't care. Some developers use VecSet(X,0) instead of
VecZeroEntries and there are many calls to VecSet(X,1.0) in PETSc but no
dedicated interface.
