[petsc-dev] TS RHS Jacobian caching (shift/scaling) in PETSc

[Stefano Zampini]

Jed,

It appears that the testsuite does not currently cover a nasty code we have
to cache RHS Jacobian computation when using the implicit interface. If I
run the testsuite with the below modification, no tests fail. (my
continuous adjoint codes instead show some reusage)

This caching mechanism is causing headaches to me and Hong when supporting
adjoints.
It would be nice if we can come up with a proper test (try to cover as much
user cases as possible) and provide a definitive fix for this (even
erroring if we don't support something, instead of silently doing the wrong
thing)

I'm pretty sure that when Jacobian computations are performed by TS itself,
using the proper matrices A,B obtained from TSGetIJacobian(ts,&A,B,...) the
caching is correct ( I have spent hours on this)

When different matrices are used, TSComputeIJacobian() without an IJacobian
callback  is buggy: note that  TSComputeIJacobian() is a public method.


diff --git a/src/ts/interface/ts.c b/src/ts/interface/ts.c
index 4b003c8f6e..02211b62ea 100644
--- a/src/ts/interface/ts.c
+++ b/src/ts/interface/ts.c
@@ -575,6 +575,7 @@ PetscErrorCode  TSComputeRHSJacobian(TS ts,PetscReal
t,Vec U,Mat A,Mat B)
   if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR ||
(ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) &&
(rhsfunction != TSComputeRHSFunctionLinear)) {
     /* restore back RHS Jacobian matrices if they have been shifted/scaled
*/
     if (A == ts->Arhs) {
+      PetscPrintf(PetscObjectComm((PetscObject)ts),"REUSE %g
%g\n",ts->rhsjacobian.shift,ts->rhsjacobian.scale);
       if (ts->rhsjacobian.shift != 0) {
         ierr = MatShift(A,-ts->rhsjacobian.shift);CHKERRQ(ierr);
       }

-- 
Stefano
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20190425/89abafec/attachment.html>