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

[Jed Brown]

Can you make a PR with a test to cover this feature?

We desperately need coverage as a CI integration.

Stefano Zampini <stefano.zampini at gmail.com> writes:

> 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