Actual source code: snestest.c
petsc-3.6.1 2015-08-06
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: }