[petsc-dev] Proposed changes to TS API

[Matthew Knepley]

On Fri, May 11, 2018 at 12:25 AM, Smith, Barry F. <bsmith at mcs.anl.gov>
wrote:

>
>
> > On May 10, 2018, at 4: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?
> >
> >> 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.
> >
> >> 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
> >
> >> 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,
>
>    What do you mean by longer critical path?
>
> > more storage required, and more data motion.
>
>   The extra storage needed is related to the size of the mass matrix
> correct? And the extra data motion is related to the size of the mass
> matrix correct?
>
>    Is the extra work coming from a needed call to MatAXPY (to combine the
> scaled mass matrix with the Jacobian) in Hong's approach? While, in theory,
> the user can avoid the MatAXPY in the current code if they have custom code
> that assembles directly the scaled mass matrix and Jacobian? But surely
> most users would not write such custom code and would themselves keep a
> copy of the mass matrix (likely constant) and use MatAXPY() to combine the
> copy of the mass matrix with the Jacobian they compute at each
> timestep/stage?  Or am I missing something?


I assemble the combined system directly.

   Matt

> And if the mass matrix is simple, won't it take a very small fraction of
> > time, thus have little gains from "optimizing it"?
>
>    Is the Function approach only theoretically much more efficient than
> Hong's approach when the mass matrix is nontrivial? That is the mass matrix
> has a nonzero structure similar to the Jacobian?
> >
> >>  *   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/20180511/95e0311b/attachment-0001.html>