Fair enough. Here you have an example. You need to checkout the branch stefano_zampini/feature-continuousadjoint and apply patch_ex23 to src/ts/examples/tutorials/ex23.c The patch just dumps the final RHS jacobian matrix of a standard TSSolve, recomputes it with by calling TSComputeRHSJacobian, and dumps it again. If you run ex23 using the RHS interface and an implicit solver, you see from the first MatView that ts->Arhs is shifted. After a call to TSComputeRHSJacobian it is ok. [szampini@localhost tutorials]$ ./ex23 -t0 0 -tf 1 -dt 0.5 -p 1 -ts_type bdf Mat Object: 1 MPI processes type: seqaij row 0: (0, 13.2008) Mat Object: 1 MPI processes type: seqaij row 0: (0, 0.7) Objective function: time [0.,1.], val 0.725152 (should be 0.724109) Vec Object: 1 MPI processes type: seq 1.45054 0.404627 -0.126697 If you apply patch_ts (which disables my fix at https://bitbucket.org/petsc/petsc/branch/stefano_zampini/feature-continuousa...), you see that the matrix is unchanged after TSComputeRHSJacobian [szampini@localhost tutorials]$ ./ex23 -t0 0 -tf 1 -dt 0.5 -p 1 -ts_type bdf Mat Object: 1 MPI processes type: seqaij row 0: (0, 13.2008) Mat Object: 1 MPI processes type: seqaij row 0: (0, 13.2008) Objective function: time [0.,1.], val 0.725152 (should be 0.724109) Vec Object: 1 MPI processes type: seq 1.45054 0.404627 -0.126697 2017-09-10 22:43 GMT+03:00 Barry Smith <[email protected]>:

Stefano,

Thanks, but I don't really want a fix. I want a code that I can run that demonstrates a problem (then we can evaluate fixes)

Barry

...

On Sep 8, 2017, at 11:51 AM, Stefano Zampini <[email protected]> wrote:

I encountered the problem while working on the adjoint code. Here is the fix I have in my branch

https://bitbucket.org/petsc/petsc/branch/stefano_zampini/ feature-continuousadjoint#Lsrc/ts/interface/ts.cT613

However, I'm not sure if this is a problem just with my code. This is why I'm asking on the list.

2017-09-08 0:02 GMT+03:00 Barry Smith <[email protected]>:

Example that demonstrates the current code is problematic in some cases?

...

On Aug 28, 2017, at 11:59 AM, Stefano Zampini < [email protected]> wrote:

I was wondering if the following code in TSComputeRHSJacobian is correct when the user pass the reuse flag via TSRHSJacobianReuse

if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) { PetscFunctionReturn(0); }

Should the above code be replaced by the following?

if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) { if (ts->rhsjacobian.reuse) { ierr = MatShift(A,-ts->rhsjacobian.shift);CHKERRQ(ierr); ierr = MatScale(A,1./ts->rhsjacobian.scale);CHKERRQ(ierr); if (B && A != B) { ierr = MatShift(B,-ts->rhsjacobian.shift);CHKERRQ(ierr); ierr = MatScale(B,1./ts->rhsjacobian.scale);CHKERRQ(ierr); } ts->rhsjacobian.shift = 0; ts->rhsjacobian.scale = 1.; } PetscFunctionReturn(0); }

-- Stefano

-- Stefano

-- Stefano