<div dir="ltr">Fair enough.<div><br></div><div>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</div><div>The patch just dumps the final RHS jacobian matrix of a standard TSSolve, recomputes it with by calling TSComputeRHSJacobian, and dumps it again.</div><div><div>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.</div><div><br></div><div>[szampini@localhost tutorials]$ ./ex23 -t0 0 -tf 1 -dt 0.5 -p 1 -ts_type bdf</div><div>Mat Object: 1 MPI processes</div><div> type: seqaij</div><div>row 0: (0, 13.2008) </div><div>Mat Object: 1 MPI processes</div><div> type: seqaij</div><div>row 0: (0, 0.7) </div><div>Objective function: time [0.,1.], val 0.725152 (should be 0.724109)</div><div>Vec Object: 1 MPI processes</div><div> type: seq</div><div>1.45054</div><div>0.404627</div><div>-0.126697</div></div><div><br></div><div><br></div><div>If you apply patch_ts (which disables my fix at <a href="https://bitbucket.org/petsc/petsc/branch/stefano_zampini/feature-continuousadjoint#Lsrc/ts/interface/ts.cT614">https://bitbucket.org/petsc/petsc/branch/stefano_zampini/feature-continuousadjoint#Lsrc/ts/interface/ts.cT614</a>), you see that the matrix is unchanged after TSComputeRHSJacobian</div><div><br></div><div><div>[szampini@localhost tutorials]$ ./ex23 -t0 0 -tf 1 -dt 0.5 -p 1 -ts_type bdf</div><div>Mat Object: 1 MPI processes</div><div> type: seqaij</div><div>row 0: (0, 13.2008) </div><div>Mat Object: 1 MPI processes</div><div> type: seqaij</div><div>row 0: (0, 13.2008) </div><div>Objective function: time [0.,1.], val 0.725152 (should be 0.724109)</div><div>Vec Object: 1 MPI processes</div><div> type: seq</div><div>1.45054</div><div>0.404627</div><div>-0.126697</div></div><div><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-09-10 22:43 GMT+03:00 Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
Stefano,<br>
<br>
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)<br>
<span class="HOEnZb"><font color="#888888"><br>
Barry<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
> On Sep 8, 2017, at 11:51 AM, Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.com</a>> wrote:<br>
><br>
> I encountered the problem while working on the adjoint code. Here is the fix I have in my branch<br>
><br>
> <a href="https://bitbucket.org/petsc/petsc/branch/stefano_zampini/feature-continuousadjoint#Lsrc/ts/interface/ts.cT613" rel="noreferrer" target="_blank">https://bitbucket.org/petsc/<wbr>petsc/branch/stefano_zampini/<wbr>feature-continuousadjoint#<wbr>Lsrc/ts/interface/ts.cT613</a><br>
><br>
> However, I'm not sure if this is a problem just with my code. This is why I'm asking on the list.<br>
><br>
><br>
><br>
> 2017-09-08 0:02 GMT+03:00 Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>>:<br>
><br>
> Example that demonstrates the current code is problematic in some cases?<br>
><br>
> > On Aug 28, 2017, at 11:59 AM, Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.com</a>> wrote:<br>
> ><br>
> ><br>
> > I was wondering if the following code in TSComputeRHSJacobian is correct when the user pass the reuse flag via TSRHSJacobianReuse<br>
> ><br>
> > if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) {<br>
> > PetscFunctionReturn(0);<br>
> > }<br>
> ><br>
> > Should the above code be replaced by the following?<br>
> ><br>
> > if (ts->rhsjacobian.time == t && (ts->problem_type == TS_LINEAR || (ts->rhsjacobian.Xid == Uid && ts->rhsjacobian.Xstate == Ustate)) && (rhsfunction != TSComputeRHSFunctionLinear)) {<br>
> > if (ts->rhsjacobian.reuse) {<br>
> > ierr = MatShift(A,-ts->rhsjacobian.<wbr>shift);CHKERRQ(ierr);<br>
> > ierr = MatScale(A,1./ts->rhsjacobian.<wbr>scale);CHKERRQ(ierr);<br>
> > if (B && A != B) {<br>
> > ierr = MatShift(B,-ts->rhsjacobian.<wbr>shift);CHKERRQ(ierr);<br>
> > ierr = MatScale(B,1./ts->rhsjacobian.<wbr>scale);CHKERRQ(ierr);<br>
> > }<br>
> > ts->rhsjacobian.shift = 0;<br>
> > ts->rhsjacobian.scale = 1.;<br>
> > }<br>
> > PetscFunctionReturn(0);<br>
> > }<br>
> ><br>
> ><br>
> > --<br>
> > Stefano<br>
><br>
><br>
><br>
><br>
> --<br>
> Stefano<br>
<br>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature">Stefano</div>
</div>