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