[petsc-users] Making a FD Jacobian

[Jed Brown]

"Mark F. Adams" <mfadams at lbl.gov> writes:

> I'm trying to make a Jacobian and have this code working:
>
>   ierr = DMCreateMatrix(dm, &A);CHKERRQ(ierr);
>
>   { /* finite difference Jacobian A */
>     SNES snes;
>     MatStructure flag;
>
>     ierr = SNESCreate(comm, &snes);CHKERRQ(ierr);
>     ierr = SNESSetDM(snes,dm);CHKERRQ(ierr);
>     ierr = SNESSetFunction(snes,NULL,ComputeFunctionLap,user);CHKERRQ(ierr);
>     ierr = SNESSetJacobian(snes,A,A,SNESComputeJacobianDefaultColor,0);CHKERRQ(ierr);
>
>     ierr = SNESComputeJacobian(snes,X,&A,&A,&flag);CHKERRQ(ierr);

The API normally requires that the residual be evaluated before the
Jacobian, so that users can compute any part of the Jacobian during
residual evaluation.  That assumption could be weakened slightly if
necessary, but if you're writing library code, that would need to be a
conscious decision.

>     /* ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
>     ierr = SNESDestroy(&snes);CHKERRQ(ierr);
>   }
>
> except I need to comment out these line in snesj2.c line ~75:
>
>   ierr = VecEqual(x1,snes->vec_sol,&solvec);CHKERRQ(ierr);

You could just skip this fast path if snes->vec_sol is NULL.

>   if (!snes->vec_rhs && solvec) {
>     ierr = MatFDColoringSetF(color,F);CHKERRQ(ierr);
>   }
>   
> The problem is that snes->vec_sol is NULL and VecEqual asserts that these are valid arguments.  How should I go about doing this?
>
> Mark
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 835 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131019/a1a6ce04/attachment.pgp>