Actual source code: snestest.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  2: #include <petsc/private/snesimpl.h>

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


 11: PetscErrorCode SNESSolve_Test(SNES snes)
 12: {
 13:   Mat            A = snes->jacobian,B;
 14:   Vec            x = snes->vec_sol,f = snes->vec_func,f1 = snes->vec_sol_update;
 16:   PetscInt       i;
 17:   PetscReal      nrm,gnorm;
 18:   SNES_Test      *neP = (SNES_Test*)snes->data;
 19:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
 20:   void           *ctx;
 21:   PetscReal      fnorm,f1norm,dnorm;

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

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

 33:   for (i=0; i<3; i++) {
 34:     void                     *functx;
 35:     static const char *const loc[] = {"user-defined state","constant state -1.0","constant state 1.0"};
 36:     PetscInt                 m,n,M,N;

 38:     if (i == 1) {
 39:       VecSet(x,-1.0);
 40:     } else if (i == 2) {
 41:       VecSet(x,1.0);
 42:     }

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

 52:     /* compute both versions of Jacobian */
 53:     SNESComputeJacobian(snes,x,A,A);

 55:     MatCreate(PetscObjectComm((PetscObject)A),&B);
 56:     MatGetSize(A,&M,&N);
 57:     MatGetLocalSize(A,&m,&n);
 58:     MatSetSizes(B,m,n,M,N);
 59:     MatSetType(B,((PetscObject)A)->type_name);
 60:     MatSetUp(B);
 61:     MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 63:     SNESGetFunction(snes,NULL,NULL,&functx);
 64:     SNESComputeJacobianDefault(snes,x,B,B,functx);
 65:     if (neP->complete_print) {
 66:       MPI_Comm    comm;
 67:       PetscViewer viewer;
 68:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Finite difference Jacobian (%s)\n",loc[i]);
 69:       PetscObjectGetComm((PetscObject)B,&comm);
 70:       PetscViewerASCIIGetStdout(comm,&viewer);
 71:       MatView(B,viewer);
 72:     }
 73:     /* compare */
 74:     MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
 75:     MatNorm(B,NORM_FROBENIUS,&nrm);
 76:     MatNorm(A,NORM_FROBENIUS,&gnorm);
 77:     if (neP->complete_print) {
 78:       MPI_Comm    comm;
 79:       PetscViewer viewer;
 80:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Hand-coded Jacobian (%s)\n",loc[i]);
 81:       PetscObjectGetComm((PetscObject)B,&comm);
 82:       PetscViewerASCIIGetStdout(comm,&viewer);
 83:       MatView(A,viewer);
 84:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Hand-coded minus finite-difference Jacobian (%s)\n",loc[i]);
 85:       MatView(B,viewer);
 86:     }
 87:     if (!gnorm) gnorm = 1; /* just in case */
 88:     PetscPrintf(PetscObjectComm((PetscObject)snes),"Norm of matrix ratio %g, difference %g (%s)\n",(double)(nrm/gnorm),(double)nrm,loc[i]);

 90:     SNESGetObjective(snes,&objective,&ctx);
 91:     if (objective) {
 92:       SNESComputeFunction(snes,x,f);
 93:       VecNorm(f,NORM_2,&fnorm);
 94:       if (neP->complete_print) {
 95:         PetscViewer viewer;
 96:         PetscPrintf(PetscObjectComm((PetscObject)snes),"Hand-coded Function (%s)\n",loc[i]);
 97:         PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
 98:         VecView(f,viewer);
 99:       }
100:       SNESObjectiveComputeFunctionDefaultFD(snes,x,f1,NULL);
101:       VecNorm(f1,NORM_2,&f1norm);
102:       if (neP->complete_print) {
103:         PetscViewer viewer;
104:         PetscPrintf(PetscObjectComm((PetscObject)snes),"Finite-difference Function (%s)\n",loc[i]);
105:         PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
106:         VecView(f1,viewer);
107:       }
108:       /* compare the two */
109:       VecAXPY(f,-1.0,f1);
110:       VecNorm(f,NORM_2,&dnorm);
111:       if (!fnorm) fnorm = 1.;
112:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Norm of function ratio %g, difference %g (%s)\n",dnorm/fnorm,dnorm,loc[i]);
113:       if (neP->complete_print) {
114:         PetscViewer viewer;
115:         PetscPrintf(PetscObjectComm((PetscObject)snes),"Difference (%s)\n",loc[i]);
116:         PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
117:         VecView(f,viewer);
118:       }
119:     }
120:     MatDestroy(&B);
121:   }

123:   /*
124:    Abort after the first iteration due to the jacobian not being valid.
125:   */

127:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESTest aborts after Jacobian test: it is NORMAL behavior.");
128:   return(0);
129: }


132: /* ------------------------------------------------------------ */
135: PetscErrorCode SNESDestroy_Test(SNES snes)
136: {

140:   PetscFree(snes->data);
141:   return(0);
142: }

146: static PetscErrorCode SNESSetFromOptions_Test(PetscOptions *PetscOptionsObject,SNES snes)
147: {
148:   SNES_Test      *ls = (SNES_Test*)snes->data;

152:   PetscOptionsHead(PetscOptionsObject,"Hand-coded Jacobian tester options");
153:   PetscOptionsBool("-snes_test_display","Display difference between hand-coded and finite difference Jacobians","None",ls->complete_print,&ls->complete_print,NULL);
154:   PetscOptionsTail();
155:   return(0);
156: }

160: PetscErrorCode SNESSetUp_Test(SNES snes)
161: {

165:   SNESSetUpMatrices(snes);
166:   return(0);
167: }

169: /* ------------------------------------------------------------ */
170: /*MC
171:       SNESTEST - Test hand-coded Jacobian against finite difference Jacobian

173:    Options Database:
174: +  -snes_type test    - use a SNES solver that evaluates the difference between hand-code and finite-difference Jacobians
175: -  -snes_test_display - display the elements of the matrix, the difference between the Jacobian approximated by finite-differencing and hand-coded Jacobian

177:    Level: intermediate

179:    Notes: This solver is not a solver and does not converge to a solution.  SNESTEST checks the Jacobian at three
180:    points: the 0, 1, and -1 solution vectors.  At each point the following is reported.

182:    Output:
183: +  difference - ||J - Jd||, the norm of the difference of the hand-coded Jacobian J and the approximate Jacobian Jd obtained by finite-differencing
184:    the residual,
185: -  ratio      - ||J - Jd||/||J||, the ratio of the norms of the above difference and the hand-coded Jacobian.

187:    Frobenius norm is used in the above throughout. After doing these three tests, it always aborts with the error message
188:    "SNESTest aborts after Jacobian test".  No other behavior is to be expected.  It may be similarly used to check if a
189:    SNES function is the gradient of an objective function set with SNESSetObjective().

191: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESUpdateCheckJacobian(), SNESNEWTONLS, SNESNEWTONTR

193: M*/
196: PETSC_EXTERN PetscErrorCode SNESCreate_Test(SNES snes)
197: {
198:   SNES_Test      *neP;

202:   snes->ops->solve          = SNESSolve_Test;
203:   snes->ops->destroy        = SNESDestroy_Test;
204:   snes->ops->setfromoptions = SNESSetFromOptions_Test;
205:   snes->ops->view           = 0;
206:   snes->ops->setup          = SNESSetUp_Test;
207:   snes->ops->reset          = 0;

209:   snes->usesksp = PETSC_FALSE;

211:   PetscNewLog(snes,&neP);
212:   snes->data          = (void*)neP;
213:   neP->complete_print = PETSC_FALSE;
214:   return(0);
215: }

219: /*@C
220:     SNESUpdateCheckJacobian - Checks each Jacobian computed by the nonlinear solver comparing the users function with a finite difference computation.

222:    Options Database:
223: +    -snes_check_jacobian - use this every time SNESSolve() is called
224: -    -snes_check_jacobian_view -  Display difference between Jacobian approximated by finite-differencing and the hand-coded Jacobian

226:    Output:
227: +  difference - ||J - Jd||, the norm of the difference of the hand-coded Jacobian J and the approximate Jacobian Jd obtained by finite-differencing
228:    the residual,
229: -  ratio      - ||J - Jd||/||J||, the ratio of the norms of the above difference and the hand-coded Jacobian.

231:    Notes:
232:    Frobenius norm is used in the above throughout.  This check is carried out every SNES iteration.

234:    Level: intermediate

236: .seealso:  SNESTEST, SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESSolve()

238: @*/
239: PetscErrorCode SNESUpdateCheckJacobian(SNES snes,PetscInt it)
240: {
241:   Mat            A = snes->jacobian,B;
242:   Vec            x = snes->vec_sol,f = snes->vec_func,f1 = snes->vec_sol_update;
244:   PetscReal      nrm,gnorm;
245:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
246:   void           *ctx;
247:   PetscReal      fnorm,f1norm,dnorm;
248:   PetscInt       m,n,M,N;
249:   PetscBool      complete_print = PETSC_FALSE;
250:   void           *functx;
251:   PetscViewer    viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));

254:   PetscOptionsHasName(((PetscObject)snes)->prefix,"-snes_check_jacobian_view",&complete_print);
255:   if (A != snes->jacobian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot check Jacobian with alternative preconditioner");

257:   PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
258:   PetscViewerASCIIPrintf(viewer,"      Testing hand-coded Jacobian, if the ratio is O(1.e-8), the hand-coded Jacobian is probably correct.\n");
259:   if (!complete_print) {
260:     PetscViewerASCIIPrintf(viewer,"      Run with -snes_check_jacobian_view [viewer][:filename][:format] to show difference of hand-coded and finite difference Jacobian.\n");
261:   }

263:   /* compute both versions of Jacobian */
264:   SNESComputeJacobian(snes,x,A,A);

266:   MatCreate(PetscObjectComm((PetscObject)A),&B);
267:   MatGetSize(A,&M,&N);
268:   MatGetLocalSize(A,&m,&n);
269:   MatSetSizes(B,m,n,M,N);
270:   MatSetType(B,((PetscObject)A)->type_name);
271:   MatSetUp(B);
272:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
273:   SNESGetFunction(snes,NULL,NULL,&functx);
274:   SNESComputeJacobianDefault(snes,x,B,B,functx);

276:   if (complete_print) {
277:     PetscViewerASCIIPrintf(viewer,"    Finite difference Jacobian\n");
278:     MatViewFromOptions(B,(PetscObject)snes,"-snes_check_jacobian_view");
279:   }
280:   /* compare */
281:   MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
282:   MatNorm(B,NORM_FROBENIUS,&nrm);
283:   MatNorm(A,NORM_FROBENIUS,&gnorm);
284:   if (complete_print) {
285:     PetscViewerASCIIPrintf(viewer,"    Hand-coded Jacobian\n");
286:     MatViewFromOptions(A,(PetscObject)snes,"-snes_check_jacobian_view");
287:     PetscViewerASCIIPrintf(viewer,"    Hand-coded minus finite difference Jacobian\n");
288:     MatViewFromOptions(B,(PetscObject)snes,"-snes_check_jacobian_view");
289:   }
290:   if (!gnorm) gnorm = 1; /* just in case */
291:   PetscViewerASCIIPrintf(viewer,"    %g = ||J - Jfd||/||J|| %g  = ||J - Jfd||\n",(double)(nrm/gnorm),(double)nrm);

293:   SNESGetObjective(snes,&objective,&ctx);
294:   if (objective) {
295:     SNESComputeFunction(snes,x,f);
296:     VecNorm(f,NORM_2,&fnorm);
297:     if (complete_print) {
298:       PetscViewerASCIIPrintf(viewer,"    Hand-coded Objective Function \n");
299:       VecView(f,viewer);
300:     }
301:     SNESObjectiveComputeFunctionDefaultFD(snes,x,f1,NULL);
302:     VecNorm(f1,NORM_2,&f1norm);
303:     if (complete_print) {
304:       PetscViewerASCIIPrintf(viewer,"    Finite-Difference Objective Function\n");
305:       VecView(f1,viewer);
306:     }
307:     /* compare the two */
308:     VecAXPY(f,-1.0,f1);
309:     VecNorm(f,NORM_2,&dnorm);
310:     if (!fnorm) fnorm = 1.;
311:     PetscViewerASCIIPrintf(viewer,"    %g = Norm of objective function ratio %g = difference\n",dnorm/fnorm,dnorm);
312:     if (complete_print) {
313:       PetscViewerASCIIPrintf(viewer,"    Difference\n");
314:       VecView(f,viewer);
315:     }
316:   }
317:   PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);

319:   MatDestroy(&B);
320:   return(0);
321: }