[petsc-dev] Proposed changes to TS API

[Matthew Knepley]

On Thu, May 10, 2018 at 5:12 PM, Jed Brown <jed at jedbrown.org> wrote:

> "Zhang, Hong" <hongzhang at anl.gov> writes:
>
> > Dear PETSc folks,
> >
> > Current TS APIs (IFunction/IJacobian+RHSFunction/RHSJacobian) were
> designed for the fully implicit formulation F(t,U,Udot) = G(t,U).
> > Shampine's paper (https://www.sciencedirect.com/science/article/pii/
> S0377042706004110?via%3Dihub<https://www.sciencedirect.com/
> science/article/pii/S0377042706004110?via=ihub>) explains some reasoning
> behind it.
> >
> > Our formulation is general enough to cover all the following common cases
> >
> >   *   Udot = G(t,U) (classic ODE)
> >   *   M Udot = G(t,U)  (ODEs/DAEs for mechanical and electronic systems)
> >   *   M(t,U) Udot = G(t,U) (PDEs)
> >
> > Yet the TS APIs provide the capability to solve both simple problems and
> complicated problems. However, we are not doing well to make TS easy to use
> and efficient especially for simple problems. Over the years, we have
> clearly seen the main drawbacks including:
> > 1. The shift parameter exposed in IJacobian is terribly confusing,
> especially to new users. Also it is not conceptually straightforward when
> using AD or finite differences on IFunction to approximate IJacobian.
>
> What isn't straightforward about AD or FD on the IFunction?  That one
> bit of chain rule?


Is it automated like SNES?


>
> > 2. It is difficult to switch from fully implicit to fully explicit.
> Users cannot use explicit methods when they provide IFunction/IJacobian.
>
> This is a real issue, but it's extremely common for PDE to have boundary
> conditions enforced as algebraic constraints, thus yielding a DAE.


Good implementations eliminate those variables :)


> > 3. The structure of mass matrix is completely invisible to TS. This
> means giving up all the opportunities to improve efficiency. For example,
> when M is constant or weekly dependent on U, we might not want to
> evaluate/update it every time IJacobian is called. If M is diagonal, the
> Jacobian can be shifted more efficiently than just using MatAXPY().
>
> I don't understand


It definitely isn't a separate entity, but I do not see that it has to be.


> > 4. Reshifting the Jacobian is unnecessarily expensive and sometimes
> buggy.
>
> Why is "reshifting" needed?  After a step is rejected and when the step
> size changes for a linear constant operator?
>
> > Consider the scenario below.
> > shift = a;
> >   TSComputeIJacobian()
> >   shift = b;
> >   TSComputeIJacobian() // with the same U and Udot as last call
> > Changing the shift parameter requires the Jacobian function to be
> evaluated again. If users provide only RHSJacobian, the Jacobian will not
> be updated/reshifted in the second call because TSComputeRHSJacobian()
> finds out that U has not been changed. This issue is fixable by adding more
> logic into the already convoluted implementation of TSComputeIJacobian(),
> but the intention here is to illustrate the cumbersomeness of current
> IJacobian and the growing complications in TSComputeIJacobian() that
> IJacobian causes.
> >
> > So I propose that we have two separate matrices dF/dUdot and dF/dU, and
> remove the shift parameter from IJacobian. dF/dU will be calculated by
> IJacobian; dF/dUdot will be calculated by a new callback function and
> default to an identity matrix if it is not provided by users. Then the
> users do not need to assemble the shifted Jacobian since TS will handle the
> shifting whenever needed. And knowing the structure of dF/dUdot (the mass
> matrix), TS will become more flexible. In particular, we will have
> >
> >   *   easy access to the unshifted Jacobian dF/dU (this simplifies the
> adjoint implementation a lot),
>
> How does this simplify the adjoint?
>
> >   *   plenty of opportunities to optimize TS when the mass matrix is
> diagonal or constant or weekly dependent on U (which accounts for almost
> all use cases in practice),
>
> But longer critical path, more storage required, and more data motion.
> And if the mass matrix is simple, won't it take a very small fraction of
> time, thus have little gains from "optimizing it"?
>

I agree here.

  Matt


> >   *   easy conversion from fully implicit to fully explicit,
> >   *   an IJacobian without shift parameter that is easy to understand
> and easy to port to other software.
>
> Note that even CVODE has an interface similar to PETSc; e.g., gamma
> parameter in CVSpilsPrecSetupFn.
>
> > Regarding the changes on the user side, most IJacobian users should not
> have problem splitting the old Jacobian if they compute dF/dUdot and dF/dU
> explicitly. If dF/dUdot is too complicated to build, matrix-free is an
> alternative option.
> >
> > While this proposal is somehow related to Barry's idea of having a
> TSSetMassMatrix() last year, I hope it provides more details for your
> information. Any of your comments and feedback would be appreciated.
> >
> > Thanks,
> > Hong
>



-- 
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/20180510/39f64f08/attachment.html>