[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