2: #include <petsc-private/snesimpl.h>

  4: typedef struct {
  5:   PetscBool  complete_print;
  6: } SNES_Test;

  8: /*
  9:      SNESSolve_Test - Tests whether a hand computed Jacobian 
 10:      matches one compute via finite differences.
 11: */
 14: PetscErrorCode SNESSolve_Test(SNES snes)
 15: {
 16:   Mat            A = snes->jacobian,B;
 17:   Vec            x = snes->vec_sol,f = snes->vec_func;
 19:   PetscInt       i;
 20:   MatStructure   flg;
 21:   PetscReal      nrm,gnorm;
 22:   SNES_Test      *neP = (SNES_Test*)snes->data;


 26:   if (A != snes->jacobian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot test with alternative preconditioner");

 28:   PetscPrintf(((PetscObject)snes)->comm,"Testing hand-coded Jacobian, if the ratio is\n");
 29:   PetscPrintf(((PetscObject)snes)->comm,"O(1.e-8), the hand-coded Jacobian is probably correct.\n");
 30:   if (!neP->complete_print) {
 31:     PetscPrintf(((PetscObject)snes)->comm,"Run with -snes_test_display to show difference\n");
 32:     PetscPrintf(((PetscObject)snes)->comm,"of hand-coded and finite difference Jacobian.\n");
 33:   }

 35:   for (i=0; i<3; i++) {
 36:     void *functx;
 37:     static const char *const loc[] = {"user-defined state","constant state -1.0","constant state 1.0"};
 38:     if (i == 1) {VecSet(x,-1.0);}
 39:     else if (i == 2) {VecSet(x,1.0);}

 41:     /* evaluate the function at this point because SNESDefaultComputeJacobianColor() assumes that the function has been evaluated and put into snes->vec_func */
 42:     SNESComputeFunction(snes,x,f);
 43:     if (snes->domainerror) {
 44:       PetscPrintf(((PetscObject)snes)->comm,"Domain error at %s\n",loc[i]);
 45:       snes->domainerror = PETSC_FALSE;
 46:       continue;
 47:     }

 49:     /* compute both versions of Jacobian */
 50:     SNESComputeJacobian(snes,x,&A,&A,&flg);
 51:     if (!i) {
 52:       PetscInt m,n,M,N;
 53:       MatCreate(((PetscObject)A)->comm,&B);
 54:       MatGetSize(A,&M,&N);
 55:       MatGetLocalSize(A,&m,&n);
 56:       MatSetSizes(B,m,n,M,N);
 57:       MatSetType(B,((PetscObject)A)->type_name);
 58:       MatSetUp(B);
 59:     }
 60:     SNESGetFunction(snes,PETSC_NULL,PETSC_NULL,&functx);
 61:     SNESDefaultComputeJacobian(snes,x,&B,&B,&flg,functx);
 62:     if (neP->complete_print) {
 63:       MPI_Comm    comm;
 64:       PetscViewer viewer;
 65:       PetscPrintf(((PetscObject)snes)->comm,"Finite difference Jacobian (%s)\n",loc[i]);
 66:       PetscObjectGetComm((PetscObject)B,&comm);
 67:       PetscViewerASCIIGetStdout(comm,&viewer);
 68:       MatView(B,viewer);
 69:     }
 70:     /* compare */
 71:     MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
 72:     MatNorm(B,NORM_FROBENIUS,&nrm);
 73:     MatNorm(A,NORM_FROBENIUS,&gnorm);
 74:     if (neP->complete_print) {
 75:       MPI_Comm    comm;
 76:       PetscViewer viewer;
 77:       PetscPrintf(((PetscObject)snes)->comm,"Hand-coded Jacobian (%s)\n",loc[i]);
 78:       PetscObjectGetComm((PetscObject)B,&comm);
 79:       PetscViewerASCIIGetStdout(comm,&viewer);
 80:       MatView(A,viewer);
 81:       PetscPrintf(((PetscObject)snes)->comm,"Hand-coded minus finite difference Jacobian (%s)\n",loc[i]);
 82:       MatView(B,viewer);
 83:     }
 84:     if (!gnorm) gnorm = 1; /* just in case */
 85:     PetscPrintf(((PetscObject)snes)->comm,"Norm of matrix ratio %g difference %g (%s)\n",(double)(nrm/gnorm),(double)nrm,loc[i]);
 86:   }
 87:   MatDestroy(&B);
 88:   /*
 89:          Return error code cause Jacobian not good
 90:   */
 91:   PetscFunctionReturn(PETSC_ERR_ARG_WRONGSTATE);
 92: }
 93: /* ------------------------------------------------------------ */
 96: PetscErrorCode SNESDestroy_Test(SNES snes)
 97: {
100:   PetscFree(snes->data);
101:   return(0);
102: }

106: static PetscErrorCode SNESSetFromOptions_Test(SNES snes)
107: {
108:   SNES_Test      *ls = (SNES_Test *)snes->data;


113:   PetscOptionsHead("Hand-coded Jacobian tester options");
114:   PetscOptionsBool("-snes_test_display","Display difference between hand-coded and finite difference Jacobians","None",ls->complete_print,&ls->complete_print,PETSC_NULL);
115:   PetscOptionsTail();
116:   return(0);
117: }

119: /* ------------------------------------------------------------ */
120: /*MC
121:       SNESTEST - Test hand-coded Jacobian against finite difference Jacobian

123:    Options Database:
124: .    -snes_test_display  Display difference between approximate and hand-coded Jacobian

126:    Level: intermediate

128: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR

130: M*/
131: EXTERN_C_BEGIN
134: PetscErrorCode  SNESCreate_Test(SNES  snes)
135: {
136:   SNES_Test      *neP;

140:   snes->ops->solve           = SNESSolve_Test;
141:   snes->ops->destroy         = SNESDestroy_Test;
142:   snes->ops->setfromoptions  = SNESSetFromOptions_Test;
143:   snes->ops->view            = 0;
144:   snes->ops->setup           = 0;
145:   snes->ops->reset           = 0;

147:   snes->usesksp             = PETSC_FALSE;

149:   PetscNewLog(snes,SNES_Test,&neP);
150:   snes->data            = (void*)neP;
151:   neP->complete_print   = PETSC_FALSE;
152:   return(0);
153: }
154: EXTERN_C_END