[petsc-users] Error when using TSSetRHSJacobian+DMSetMatrixPreallocateOnly

[Vijay S. Mahadevan]

All,

When I call DMSetMatrixPreallocateOnly and then TSSetRHSJacobian
without a matrix input argument (NULL instead), then the code crashes
during TSSolve and more specifically during MatDuplicate. This is
replicable with src/ts/examples/tutorials/ex2.c. Here's the error with
stack trace.

[0]PETSC ERROR: --------------------- Error Message
------------------------------------
[0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Not for unassembled matrix!
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: Petsc Development GIT revision:
0799ffdeffb47ef0458c1305293fae28c4d2cd92  GIT Date: 2013-08-01
06:36:16 +0800
[0]PETSC ERROR: See docs/changes/index.html for recent updates.
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
[0]PETSC ERROR: See docs/index.html for manual pages.
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: ./ex2 on a arch-darwin-cxx-debug named Tachyon.local
by mahadevan Tue Sep 17 21:26:54 2013
[0]PETSC ERROR: Libraries linked from /opt/petsc-dbg/lib
[0]PETSC ERROR: Configure run at Wed Jul 31 18:03:46 2013
[0]PETSC ERROR: Configure options --prefix=/opt/petsc-dbg
PETSC_ARCH=arch-darwin-cxx-debug --download-scalapack=1
--download-blacs=1
--with-mumps-dir=/Users/mahadevan/source/MUMPS_4.10.0-p3
--download-hypre=1 --with-metis=1 --with-parmetis=1
--known-mpi-shared-libraries=1
--with-blas-lapack-lib=/System/Library/Frameworks/vecLib.framework/vecLib
--with-c++-support=1 --with-c-support=1 --with-cc=mpicc
--with-clanguage=C++ --with-moab-dir=/opt/moab
--with-dynamic-loading=1 --with-fc=mpif90 --with-fortran=1
--with-mpi=1 --with-shared-libraries=1 --with-valgrind=1
--with-valgrind-dir=/opt/local --with-cc=mpicc --with-cxx=mpicxx
COPTFLAGS="-g -fPIC" CXXOPTFLAGS="-g -fPIC" FOPTFLAGS="-g -fPIC"
--with-metis-dir=/usr/local --with-parmetis-dir=/usr/local
--with-netcdf-dir=/usr/local --with-zoltan=1
--with-zoltan-lib="-L/usr/local/lib -lptscotch -lscotch -lscotchmetis
-lscotcherr -lscotcherrexit -lzoltan"
--with-zoltan-include=/usr/local/include --with-hdf5
--with-hdf5-include=/usr/local/include
--with-hdf5-lib="-L/usr/local/lib -lhdf5_fortran -lhdf5"
--with-netcdf-dir=/usr/local --with-cmake=/opt/local/bin/cmake
--with-x-dir=/opt/local
[0]PETSC ERROR:
------------------------------------------------------------------------
[0]PETSC ERROR: MatDuplicate() line 4109 in
/Users/mahadevan/source/petsc-dev/src/mat/interface/matrix.c
[0]PETSC ERROR: TSGetRHSMats_Private() line 568 in
/Users/mahadevan/source/petsc-dev/src/ts/interface/ts.c
[0]PETSC ERROR: TSComputeIJacobian() line 740 in
/Users/mahadevan/source/petsc-dev/src/ts/interface/ts.c
[0]PETSC ERROR: SNESTSFormJacobian_Theta() line 307 in
/Users/mahadevan/source/petsc-dev/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: SNESTSFormJacobian() line 3446 in
/Users/mahadevan/source/petsc-dev/src/ts/interface/ts.c
[0]PETSC ERROR: SNESComputeJacobian() line 2227 in
/Users/mahadevan/source/petsc-dev/src/snes/interface/snes.c
[0]PETSC ERROR: SNESSolve_NEWTONLS() line 233 in
/Users/mahadevan/source/petsc-dev/src/snes/impls/ls/ls.c
[0]PETSC ERROR: SNESSolve() line 3789 in
/Users/mahadevan/source/petsc-dev/src/snes/interface/snes.c
[0]PETSC ERROR: TSStep_Theta() line 183 in
/Users/mahadevan/source/petsc-dev/src/ts/impls/implicit/theta/theta.c
[0]PETSC ERROR: TSStep() line 2507 in
/Users/mahadevan/source/petsc-dev/src/ts/interface/ts.c
[0]PETSC ERROR: TSSolve() line 2632 in
/Users/mahadevan/source/petsc-dev/src/ts/interface/ts.c
[0]PETSC ERROR: main() line 196 in
/Users/mahadevan/source/petsc-dev/src/ts/examples/tutorials/ex2.c

I hope this is not supposed to happen. The fix seems to be to call
TSSetRHSJacobian with a user created matrix instead of letting DM
create one, which rectifies this issue. This behavior feels
inconsistent and is there a way to fix this.

Vijay

PS: Just replace the following lines in ex2.c to replicate the error.

/*  ierr = TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);CHKERRQ(ierr); */
  ierr = TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr);
  ierr = DMSetMatrixPreallocateOnly(appctx.da,PETSC_TRUE);CHKERRQ(ierr);