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