[petsc-users] TSDefaultComputeJacobian

[Barry Smith]

On Jul 21, 2010, at 8:43 PM, Xuan YU wrote:

> 
> 
> Begin forwarded message:
> 
>> From: Jed Brown <jed at 59A2.org>
>> Date: July 21, 2010 8:04:17 PM EDT
>> To: Xuan YU <xxy113 at psu.edu>, PETSc users list <petsc-users at mcs.anl.gov>
>> Subject: Re: [petsc-users] TSDefaultComputeJacobian
>> 
>> On Wed, 21 Jul 2010 19:52:23 -0400, Xuan YU <xxy113 at psu.edu> wrote:
>>> Hi,
>>> 
>>> In SNES, Petsc can compute Jacobian by:
>>> 
>>> ierr =  
>>> SNESSetJacobian 
>>> (snes,J,J,SNESDefaultComputeJacobian,PETSC_NULL);CHKERRQ(ierr);
>>> 
>>> How about TS? I don't want to provide Jacobian function by myself, but  
>>> I didn't find something like TSDefaultComputeJacobian.
>> 
>> You can use -snes_fd (or -snes_mf or -snes_mf_operator) with TS.
>> Programmatically,
>> 
>>  TSGetSNES(ts,&snes);
>>  SNESSetJacobian(snes,A,B,SNESDefaultComputeJacobian,ts);
>> 
> 
> I have already tried this function.
> 
> I want to use Finite Difference Jacobian Approximations
> 
> So, before I call the function: MatGetColoring(J,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr); I should assemble J. I tried to use SNESComputeJacobian(snes,x,&J,&J,&flag) to assemble J.
> 
> Here is the code
> 
> 
> 
>         ierr = MatCreateSeqAIJ(PETSC_COMM_SELF,N,N,10,PETSC_NULL,&J);CHKERRQ(ierr);
>         ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
>         ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
>         ierr = TSMonitorSet(ts,Monitor,&appctx,PETSC_NULL);CHKERRQ(ierr);
>         ierr = TSSetSolution(ts,CV_Y);CHKERRQ(ierr);
>         ierr = TSSetRHSFunction(ts,f,&appctx);CHKERRQ(ierr);
>         ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
>         ierr = TSGetSNES(ts,&ts_snes);
>         ierr= SNESSetJacobian(ts_snes,J,J,SNESDefaultComputeJacobian,PETSC_NULL);CHKERRQ(ierr);
>         MatStructure  flag;
>         ierr = SNESComputeJacobian(ts_snes,CV_Y,&J,&J,&flag);CHKERRQ(ierr);

   Don't call this here; SNES is not ready to compute Jacobians so cannot. Just call SNESDefaultComputeJacobian() with appropriate arguments directly.

>         ierr = MatGetColoring(J,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr);

   This will not do what you want here. The default Jacobian is dense so the coloring is one color per column so the coloring does not help at all.

   Just bite the bullet and figure out the nonzero pattern of the Jacobian yourself and don't try to use SNESDefaultComputeJacobian() to "cheat" and figure it out for you. If you know the code in your RHS function then you know the nonzero pattern of the Jacobian. It is just a matter of writing the code that produces it.

   Barry

>         ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr);
>         ISColoringDestroy(iscoloring);
> 
>         ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))f,(void*)&appctx);CHKERRQ(ierr);
>         ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr);
>         ierr = TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
> 
> But I got the error:
> [0]PETSC ERROR: --------------------- Error Message ------------------------------------
> [0]PETSC ERROR: Object is in wrong state!
> [0]PETSC ERROR: Must call SNESSetFunction() before SNESComputeFunction(), likely called from SNESSolve().!
> [0]PETSC ERROR: ------------------------------------------------------------------------
> [0]PETSC ERROR: Petsc Development HG revision: fdc7be12de37b8a400b6e561de9850bcc7e79f4f  HG Date: Fri Jul 02 21:17:36 2010 -0500
> [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: ./pihm on a arch-linu named kernighan.aset.psu.edu by xxy113 Wed Jul 21 21:35:01 2010
> [0]PETSC ERROR: Libraries linked from /gpfs/home/xxy113/soft/petsc-dev/arch-linux-gnu-c-debug/lib
> [0]PETSC ERROR: Configure run at Sun Jul  4 14:09:16 2010
> [0]PETSC ERROR: Configure options --download-f-blas-lapack=1 --download-mpich=1
> [0]PETSC ERROR: ------------------------------------------------------------------------
> [0]PETSC ERROR: SNESComputeFunction() line 1180 in src/snes/interface/snes.c
> [0]PETSC ERROR: SNESDefaultComputeJacobian() line 77 in src/snes/interface/snesj.c
> [0]PETSC ERROR: SNESComputeJacobian() line 1264 in src/snes/interface/snes.c
> [0]PETSC ERROR: VecRestoreArrayPrivate3() line 234 in /work/petsc/pihm.c
> 
> 
> 
>> Jed
>> 
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100721/bc56e2d2/attachment-0001.htm>