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

Stefano Zampini stefano.zampini at gmail.com
Thu Apr 25 14:08:47 CDT 2019


I can try to excerpt a test case from PETScOpt an add it to PETSc.
On the other hand, I really don't know what are all the possible use cases
we need to cover, and if we really want to maintain the current caching
infrastructure.
I would be inclined to throw it away (I should still take a look at how
many times I am actually reusing the cached matrices)

It's unbelievable how many things are not under coverage. On average I
think I personally fix 10 bugs per week in PETSc.
For sure, the 5-point stencil Laplacian is bugfree ;-)

Using petsc4py could be an easy way to extend coverage very quickly.

Il giorno gio 25 apr 2019 alle ore 18:06 Jed Brown <jed at jedbrown.org> ha
scritto:

> 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
>


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


More information about the petsc-dev mailing list