[petsc-dev] TSComputeRHSJacobian

[Stefano Zampini]

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 at 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-continuousadjoint#Lsrc/ts/interface/ts.cT614),
you see that the matrix is unchanged after TSComputeRHSJacobian

[szampini at 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 <bsmith at mcs.anl.gov>:

>
>   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 <stefano.zampini at gmail.com>
> 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 <bsmith at mcs.anl.gov>:
> >
> >   Example that demonstrates the current code is problematic in some
> cases?
> >
> > > On Aug 28, 2017, at 11:59 AM, Stefano Zampini <
> stefano.zampini at gmail.com> 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170914/2030287b/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: patch_ex23
Type: application/octet-stream
Size: 1169 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170914/2030287b/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: patch_ts
Type: application/octet-stream
Size: 678 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20170914/2030287b/attachment-0001.obj>