[petsc-dev] -snes_mf_operator fails with TS
Barry Smith
bsmith at mcs.anl.gov
Mon Nov 11 18:18:58 CST 2013
I took at look at possibly fixing this but the following code (and the routines it calls) are so delicately balanced I had no hope. It seems to me a new factorization has to be possible to eliminate the complexity for something that really isn’t all that complicated.
Barry
PetscErrorCode TSComputeIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal shift,Mat *A,Mat *B,MatStructure *flg,PetscBool imex)
{
PetscErrorCode ierr;
TSIJacobian ijacobian;
TSRHSJacobian rhsjacobian;
DM dm;
void *ctx;
PetscFunctionBegin;
PetscValidHeaderSpecific(ts,TS_CLASSID,1);
PetscValidHeaderSpecific(U,VEC_CLASSID,3);
PetscValidHeaderSpecific(Udot,VEC_CLASSID,4);
PetscValidPointer(A,6);
PetscValidHeaderSpecific(*A,MAT_CLASSID,6);
PetscValidPointer(B,7);
PetscValidHeaderSpecific(*B,MAT_CLASSID,7);
PetscValidPointer(flg,8);
ierr = TSGetDM(ts,&dm);CHKERRQ(ierr);
ierr = DMTSGetIJacobian(dm,&ijacobian,&ctx);CHKERRQ(ierr);
ierr = DMTSGetRHSJacobian(dm,&rhsjacobian,NULL);CHKERRQ(ierr);
if (!rhsjacobian && !ijacobian) SETERRQ(PetscObjectComm((PetscObject)ts),PETSC_ERR_USER,"Must call TSSetRHSJacobian() and / or TSSetIJacobian()");
*flg = SAME_NONZERO_PATTERN; /* In case we're solving a linear problem in which case it wouldn't get initialized below. */
ierr = PetscLogEventBegin(TS_JacobianEval,ts,U,*A,*B);CHKERRQ(ierr);
if (ijacobian) {
*flg = DIFFERENT_NONZERO_PATTERN;
PetscStackPush("TS user implicit Jacobian");
ierr = (*ijacobian)(ts,t,U,Udot,shift,A,B,flg,ctx);CHKERRQ(ierr);
PetscStackPop;
/* make sure user returned a correct Jacobian and preconditioner */
PetscValidHeaderSpecific(*A,MAT_CLASSID,4);
PetscValidHeaderSpecific(*B,MAT_CLASSID,5);
}
if (imex) {
if (!ijacobian) { /* system was written as Udot = G(t,U) */
ierr = MatZeroEntries(*A);CHKERRQ(ierr);
ierr = MatShift(*A,shift);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatZeroEntries(*B);CHKERRQ(ierr);
ierr = MatShift(*B,shift);CHKERRQ(ierr);
}
*flg = SAME_PRECONDITIONER;
}
} else {
Mat Arhs = NULL,Brhs = NULL;
MatStructure flg2;
if (rhsjacobian) {
ierr = TSGetRHSMats_Private(ts,&Arhs,&Brhs);CHKERRQ(ierr);
ierr = TSComputeRHSJacobian(ts,t,U,&Arhs,&Brhs,&flg2);CHKERRQ(ierr);
}
if (Arhs == *A) { /* No IJacobian, so we only have the RHS matrix */
ts->rhsjacobian.scale = -1;
ts->rhsjacobian.shift = shift;
ierr = MatScale(*A,-1);CHKERRQ(ierr);
ierr = MatShift(*A,shift);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatScale(*B,-1);CHKERRQ(ierr);
ierr = MatShift(*B,shift);CHKERRQ(ierr);
}
} else if (Arhs) { /* Both IJacobian and RHSJacobian */
MatStructure axpy = DIFFERENT_NONZERO_PATTERN;
if (!ijacobian) { /* No IJacobian provided, but we have a separate RHS matrix */
ierr = MatZeroEntries(*A);CHKERRQ(ierr);
ierr = MatShift(*A,shift);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatZeroEntries(*B);CHKERRQ(ierr);
ierr = MatShift(*B,shift);CHKERRQ(ierr);
}
}
ierr = MatAXPY(*A,-1,Arhs,axpy);CHKERRQ(ierr);
if (*A != *B) {
ierr = MatAXPY(*B,-1,Brhs,axpy);CHKERRQ(ierr);
}
*flg = PetscMin(*flg,flg2);
}
}
On Nov 11, 2013, at 10:56 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> It is extremely frustrating that something like -snes_mf_operator which has been around since the beginning of time doesn’t work with TS
>
> Barry
>
> On Nov 11, 2013, at 10:46 AM, Jay J. Billings <billingsjj at ornl.gov> wrote:
>
>> Barry,
>>
>> Adding the -snes_mf_operator option fails with the following output
>>
>> 0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------
>> [0]PETSC ERROR: No support for this operation for this object type!
>> [0]PETSC ERROR: Not written for this matrix type!
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Petsc Development GIT revision:
>> a965ca046084fa53248a41da989a0a62cb6266ea GIT Date: 2013-11-10 14:25:12
>> -0600
>> [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: ./xolotl on a arch-linux2-c-opt named
>> antecessor.ornl.gov by bkj Mon Nov 11 11:45:31 2013
>> [0]PETSC ERROR: Libraries linked from /opt/petsc-latest_mpich-3.0.1/lib
>> [0]PETSC ERROR: Configure run at Mon Nov 11 09:36:40 2013
>> [0]PETSC ERROR: Configure options --prefix=/opt/petsc-latest_mpich-3.0.1
>> --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif77 --with-debugging=no
>> --download-f-blas-lapack=1 --FOPTFLAGS= --with-shared-libraries=1
>> --download-hypre=yes --with-debugging=0
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: MatDuplicate() line 4125 in src/mat/interface/matrix.c
>> [0]PETSC ERROR: TSGetRHSMats_Private() line 583 in src/ts/interface/ts.c
>> [0]PETSC ERROR: TSComputeIJacobian() line 760 in src/ts/interface/ts.c
>> [0]PETSC ERROR: SNESTSFormJacobian_ARKIMEX() line 1057 in
>> src/ts/impls/arkimex/arkimex.c
>> [0]PETSC ERROR: SNESTSFormJacobian() line 3543 in src/ts/interface/ts.c
>> [0]PETSC ERROR: SNESComputeJacobian() line 2241 in src/snes/interface/snes.c
>> [0]PETSC ERROR: SNESSolve_NEWTONLS() line 231 in src/snes/impls/ls/ls.c
>> [0]PETSC ERROR: SNESSolve() line 3809 in src/snes/interface/snes.c
>> [0]PETSC ERROR: TSStep_ARKIMEX() line 777 in src/ts/impls/arkimex/arkimex.c
>> [0]PETSC ERROR: TSStep() line 2603 in src/ts/interface/ts.c
>> [0]PETSC ERROR: TSSolve() line 2728 in src/ts/interface/ts.c
>> [0]PETSC ERROR: checkPetscError() line 61 in
>> /home/bkj/research/xolotl/xolotl_workspace/xolotl-Source at xolotl/xolotlSolver/PetscSolver.cpp
>>
>>
>> Jay
>>
>> ------------------------------------------------------------------------------
>> November Webinars for C, C++, Fortran Developers
>> Accelerate application performance with scalable programming models. Explore
>> techniques for threading, error checking, porting, and tuning. Get the most
>> from the latest Intel processors and coprocessors. See abstracts and register
>> http://pubads.g.doubleclick.net/gampad/clk?id=60136231&iu=/4140/ostg.clktrk
>> _______________________________________________
>> Xolotl-psi-development mailing list
>> Xolotl-psi-development at lists.sourceforge.net
>> https://lists.sourceforge.net/lists/listinfo/xolotl-psi-development
>
More information about the petsc-dev
mailing list