<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
</head>
<body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">
<div class="">Dear PETSc folks,</div>
<div class=""><br class="">
</div>
<div class="">Current TS APIs (IFunction/IJacobian+RHSFunction/RHSJacobian) were designed for the fully implicit formulation F(t,U,Udot) = G(t,U).</div>
<div class="">Shampine's paper (<a href="https://www.sciencedirect.com/science/article/pii/S0377042706004110?via=ihub" class="">https://www.sciencedirect.com/science/article/pii/S0377042706004110?via%3Dihub</a>) explains some reasoning behind it.</div>
<div class=""><br class="">
</div>
<div class="">Our formulation is general enough to cover all the following common cases</div>
<div class="">
<ul class="MailOutline">
<li class="">Udot = G(t,U) (classic ODE)</li><li class="">M Udot = G(t,U) (ODEs/DAEs for mechanical and electronic systems)</li><li class="">M(t,U) Udot = G(t,U) (PDEs)</li></ul>
</div>
<span class="">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:<br class="">
</span>
<div class=""><span class="">
<div class=""><span class="Apple-tab-span" style="white-space:pre"></span>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.</div>
<div class=""><span class="Apple-tab-span" style="white-space:pre"></span>2. It is difficult to switch from fully implicit to fully explicit. Users cannot use explicit methods when they provide IFunction/IJacobian.</div>
<div class=""><span class="Apple-tab-span" style="white-space:pre"></span>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().<br class="">
</div>
<div class=""><span class="Apple-tab-span" style="white-space:pre"></span>4. Reshifting the Jacobian is unnecessarily expensive and sometimes buggy.</div>
<div class="">Consider the scenario below.</div>
<div class=""><span class="Apple-tab-span" style="white-space:pre"></span>shift = a;</div>
</span></div>
<div class=""> <span class="Apple-tab-span" style="white-space:pre"></span>TSComputeIJacobian()</div>
<div class=""> <span class="Apple-tab-span" style="white-space:pre"></span>shift = b;</div>
<div class=""> <span class="Apple-tab-span" style="white-space:pre"></span>TSComputeIJacobian() // with the same U and Udot as last call</div>
<div class="">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.</div>
<div class=""><br class="">
</div>
<div class="">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</div>
<div class="">
<ul class="MailOutline">
<li class="">easy access to the unshifted Jacobian dF/dU (this simplifies the adjoint implementation a lot),</li><li class="">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),</li><li class="">easy conversion from fully implicit to fully explicit,</li><li class="">an IJacobian without shift parameter that is easy to understand and easy to port to other software.</li></ul>
<div class=""><br class="">
</div>
</div>
<div class="">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.</div>
<div class=""><br class="">
</div>
<div class="">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.</div>
<div class=""><br class="">
</div>
<div class="">Thanks,</div>
<div class="">Hong</div>
</body>
</html>