Actual source code: ex3.c
2: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
3: This example employs a user-defined monitoring routine and optionally a user-defined\n\
4: routine to check candidate iterates produced by line search routines. This code also\n\
5: demonstrates use of the macro __FUNCT__ to define routine names for use in error handling.\n\
6: The command line options include:\n\
7: -check_iterates : activate checking of iterates\n\
8: -check_tol <tol>: set tolerance for iterate checking\n\n";
10: /*T
11: Concepts: SNES^basic parallel example
12: Concepts: SNES^setting a user-defined monitoring routine
13: Concepts: error handling^using the macro __FUNCT__ to define routine names;
14: Processors: n
15: T*/
17: /*
18: Include "petscdraw.h" so that we can use distributed arrays (DAs).
19: Include "petscdraw.h" so that we can use PETSc drawing routines.
20: Include "petscsnes.h" so that we can use SNES solvers. Note that this
21: file automatically includes:
22: petsc.h - base PETSc routines petscvec.h - vectors
23: petscsys.h - system routines petscmat.h - matrices
24: petscis.h - index sets petscksp.h - Krylov subspace methods
25: petscviewer.h - viewers petscpc.h - preconditioners
26: petscksp.h - linear solvers
27: */
29: #include petscda.h
30: #include petscsnes.h
32: /*
33: User-defined routines. Note that immediately before each routine below,
34: we define the macro __FUNCT__ to be a string containing the routine name.
35: If defined, this macro is used in the PETSc error handlers to provide a
36: complete traceback of routine names. All PETSc library routines use this
37: macro, and users can optionally employ it as well in their application
38: codes. Note that users can get a traceback of PETSc errors regardless of
39: whether they define __FUNCT__ in application codes; this macro merely
40: provides the added traceback detail of the application routine names.
41: */
42: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
43: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
44: PetscErrorCode FormInitialGuess(Vec);
45: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void *);
46: PetscErrorCode PreCheck(SNES,Vec,Vec,void*,PetscTruth *);
47: PetscErrorCode PostCheck(SNES,Vec,Vec,Vec,void*,PetscTruth *,PetscTruth *);
49: /*
50: User-defined application context
51: */
52: typedef struct {
53: DA da; /* distributed array */
54: Vec F; /* right-hand-side of PDE */
55: PetscMPIInt rank; /* rank of processor */
56: PetscMPIInt size; /* size of communicator */
57: PetscReal h; /* mesh spacing */
58: } ApplicationCtx;
60: /*
61: User-defined context for monitoring
62: */
63: typedef struct {
64: PetscViewer viewer;
65: } MonitorCtx;
67: /*
68: User-defined context for checking candidate iterates that are
69: determined by line search methods
70: */
71: typedef struct {
72: Vec last_step; /* previous iterate */
73: PetscReal tolerance; /* tolerance for changes between successive iterates */
74: ApplicationCtx *user;
75: } StepCheckCtx;
79: int main(int argc,char **argv)
80: {
81: SNES snes; /* SNES context */
82: Mat J; /* Jacobian matrix */
83: ApplicationCtx ctx; /* user-defined context */
84: Vec x,r,U,F; /* vectors */
85: MonitorCtx monP; /* monitoring context */
86: StepCheckCtx checkP; /* step-checking context */
87: PetscTruth step_check; /* flag indicating whether we're checking
88: candidate iterates */
89: PetscScalar xp,*FF,*UU,none = -1.0;
91: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
92: PetscReal abstol,rtol,stol,norm;
94: PetscInitialize(&argc,&argv,(char *)0,help);
95: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
96: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
97: PetscOptionsGetInt(PETSC_NULL,"-n",&N,PETSC_NULL);
98: ctx.h = 1.0/(N-1);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create nonlinear solver context
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: SNESCreate(PETSC_COMM_WORLD,&snes);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create vector data structures; set function evaluation routine
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Create distributed array (DA) to manage parallel grid and vectors
112: */
113: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,N,1,1,PETSC_NULL,&ctx.da);
115: /*
116: Extract global and local vectors from DA; then duplicate for remaining
117: vectors that are the same types
118: */
119: DACreateGlobalVector(ctx.da,&x);
120: VecDuplicate(x,&r);
121: VecDuplicate(x,&F); ctx.F = F;
122: VecDuplicate(x,&U);
124: /*
125: Set function evaluation routine and vector. Whenever the nonlinear
126: solver needs to compute the nonlinear function, it will call this
127: routine.
128: - Note that the final routine argument is the user-defined
129: context that provides application-specific data for the
130: function evaluation routine.
131: */
132: SNESSetFunction(snes,r,FormFunction,&ctx);
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Create matrix data structure; set Jacobian evaluation routine
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: MatCreate(PETSC_COMM_WORLD,&J);
139: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
140: MatSetFromOptions(J);
142: /*
143: Set Jacobian matrix data structure and default Jacobian evaluation
144: routine. Whenever the nonlinear solver needs to compute the
145: Jacobian matrix, it will call this routine.
146: - Note that the final routine argument is the user-defined
147: context that provides application-specific data for the
148: Jacobian evaluation routine.
149: */
150: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Customize nonlinear solver; set runtime options
154: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
156: /*
157: Set an optional user-defined monitoring routine
158: */
159: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
160: SNESMonitorSet(snes,Monitor,&monP,0);
162: /*
163: Set names for some vectors to facilitate monitoring (optional)
164: */
165: PetscObjectSetName((PetscObject)x,"Approximate Solution");
166: PetscObjectSetName((PetscObject)U,"Exact Solution");
168: /*
169: Set SNES/KSP/KSP/PC runtime options, e.g.,
170: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
171: */
172: SNESSetFromOptions(snes);
174: /*
175: Set an optional user-defined routine to check the validity of candidate
176: iterates that are determined by line search methods
177: */
178: PetscOptionsHasName(PETSC_NULL,"-post_check_iterates",&step_check);
179: if (step_check) {
180: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
181: SNESLineSearchSetPostCheck(snes,PostCheck,&checkP);
182: VecDuplicate(x,&(checkP.last_step));
183: checkP.tolerance = 1.0;
184: checkP.user = &ctx;
185: PetscOptionsGetReal(PETSC_NULL,"-check_tol",&checkP.tolerance,PETSC_NULL);
186: }
188: PetscOptionsHasName(PETSC_NULL,"-pre_check_iterates",&step_check);
189: if (step_check) {
190: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
191: SNESLineSearchSetPreCheck(snes,PreCheck,&checkP);
192: }
194: /*
195: Print parameters used for convergence testing (optional) ... just
196: to demonstrate this routine; this information is also printed with
197: the option -snes_view
198: */
199: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
200: PetscPrintf(PETSC_COMM_WORLD,"atol=%G, rtol=%G, stol=%G, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Initialize application:
204: Store right-hand-side of PDE and exact solution
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207: /*
208: Get local grid boundaries (for 1-dimensional DA):
209: xs, xm - starting grid index, width of local grid (no ghost points)
210: */
211: DAGetCorners(ctx.da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
213: /*
214: Get pointers to vector data
215: */
216: DAVecGetArray(ctx.da,F,&FF);
217: DAVecGetArray(ctx.da,U,&UU);
219: /*
220: Compute local vector entries
221: */
222: xp = ctx.h*xs;
223: for (i=xs; i<xs+xm; i++) {
224: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
225: UU[i] = xp*xp*xp;
226: xp += ctx.h;
227: }
229: /*
230: Restore vectors
231: */
232: DAVecRestoreArray(ctx.da,F,&FF);
233: DAVecRestoreArray(ctx.da,U,&UU);
235: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: Evaluate initial guess; then solve nonlinear system
237: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239: /*
240: Note: The user should initialize the vector, x, with the initial guess
241: for the nonlinear solver prior to calling SNESSolve(). In particular,
242: to employ an initial guess of zero, the user should explicitly set
243: this vector to zero by calling VecSet().
244: */
245: FormInitialGuess(x);
246: SNESSolve(snes,PETSC_NULL,x);
247: SNESGetIterationNumber(snes,&its);
248: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n\n",its);
250: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: Check solution and clean up
252: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254: /*
255: Check the error
256: */
257: VecAXPY(x,none,U);
258: VecNorm(x,NORM_2,&norm);
259: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %D\n",norm,its);
261: /*
262: Free work space. All PETSc objects should be destroyed when they
263: are no longer needed.
264: */
265: PetscViewerDestroy(monP.viewer);
266: if (step_check) {VecDestroy(checkP.last_step);}
267: VecDestroy(x);
268: VecDestroy(r);
269: VecDestroy(U);
270: VecDestroy(F);
271: MatDestroy(J);
272: SNESDestroy(snes);
273: DADestroy(ctx.da);
274: PetscFinalize();
276: return(0);
277: }
278: /* ------------------------------------------------------------------- */
281: /*
282: FormInitialGuess - Computes initial guess.
284: Input/Output Parameter:
285: . x - the solution vector
286: */
287: PetscErrorCode FormInitialGuess(Vec x)
288: {
290: PetscScalar pfive = .50;
293: VecSet(x,pfive);
294: return(0);
295: }
296: /* ------------------------------------------------------------------- */
299: /*
300: FormFunction - Evaluates nonlinear function, F(x).
302: Input Parameters:
303: . snes - the SNES context
304: . x - input vector
305: . ctx - optional user-defined context, as set by SNESSetFunction()
307: Output Parameter:
308: . f - function vector
310: Note:
311: The user-defined context can contain any application-specific
312: data needed for the function evaluation.
313: */
314: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
315: {
316: ApplicationCtx *user = (ApplicationCtx*) ctx;
317: DA da = user->da;
318: PetscScalar *xx,*ff,*FF,d;
320: PetscInt i,M,xs,xm;
321: Vec xlocal;
324: DAGetLocalVector(da,&xlocal);
325: /*
326: Scatter ghost points to local vector, using the 2-step process
327: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
328: By placing code between these two statements, computations can
329: be done while messages are in transition.
330: */
331: DAGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
332: DAGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
334: /*
335: Get pointers to vector data.
336: - The vector xlocal includes ghost point; the vectors x and f do
337: NOT include ghost points.
338: - Using DAVecGetArray() allows accessing the values using global ordering
339: */
340: DAVecGetArray(da,xlocal,&xx);
341: DAVecGetArray(da,f,&ff);
342: DAVecGetArray(da,user->F,&FF);
344: /*
345: Get local grid boundaries (for 1-dimensional DA):
346: xs, xm - starting grid index, width of local grid (no ghost points)
347: */
348: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
349: DAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
350: PETSC_NULL,PETSC_NULL,PETSC_NULL);
352: /*
353: Set function values for boundary points; define local interior grid point range:
354: xsi - starting interior grid index
355: xei - ending interior grid index
356: */
357: if (xs == 0) { /* left boundary */
358: ff[0] = xx[0];
359: xs++;xm--;
360: }
361: if (xs+xm == M) { /* right boundary */
362: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
363: xm--;
364: }
366: /*
367: Compute function over locally owned part of the grid (interior points only)
368: */
369: d = 1.0/(user->h*user->h);
370: for (i=xs; i<xs+xm; i++) {
371: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
372: }
374: /*
375: Restore vectors
376: */
377: DAVecRestoreArray(da,xlocal,&xx);
378: DAVecRestoreArray(da,f,&ff);
379: DAVecRestoreArray(da,user->F,&FF);
380: DARestoreLocalVector(da,&xlocal);
381: return(0);
382: }
383: /* ------------------------------------------------------------------- */
386: /*
387: FormJacobian - Evaluates Jacobian matrix.
389: Input Parameters:
390: . snes - the SNES context
391: . x - input vector
392: . dummy - optional user-defined context (not used here)
394: Output Parameters:
395: . jac - Jacobian matrix
396: . B - optionally different preconditioning matrix
397: . flag - flag indicating matrix structure
398: */
399: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *ctx)
400: {
401: ApplicationCtx *user = (ApplicationCtx*) ctx;
402: PetscScalar *xx,d,A[3];
404: PetscInt i,j[3],M,xs,xm;
405: DA da = user->da;
408: /*
409: Get pointer to vector data
410: */
411: DAVecGetArray(da,x,&xx);
412: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
414: /*
415: Get range of locally owned matrix
416: */
417: DAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
418: PETSC_NULL,PETSC_NULL,PETSC_NULL);
420: /*
421: Determine starting and ending local indices for interior grid points.
422: Set Jacobian entries for boundary points.
423: */
425: if (xs == 0) { /* left boundary */
426: i = 0; A[0] = 1.0;
427: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
428: xs++;xm--;
429: }
430: if (xs+xm == M) { /* right boundary */
431: i = M-1; A[0] = 1.0;
432: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
433: xm--;
434: }
436: /*
437: Interior grid points
438: - Note that in this case we set all elements for a particular
439: row at once.
440: */
441: d = 1.0/(user->h*user->h);
442: for (i=xs; i<xs+xm; i++) {
443: j[0] = i - 1; j[1] = i; j[2] = i + 1;
444: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
445: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
446: }
448: /*
449: Assemble matrix, using the 2-step process:
450: MatAssemblyBegin(), MatAssemblyEnd().
451: By placing code between these two statements, computations can be
452: done while messages are in transition.
454: Also, restore vector.
455: */
457: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
458: DAVecRestoreArray(da,x,&xx);
459: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
461: *flag = SAME_NONZERO_PATTERN;
462: return(0);
463: }
464: /* ------------------------------------------------------------------- */
467: /*
468: Monitor - Optional user-defined monitoring routine that views the
469: current iterate with an x-window plot. Set by SNESMonitorSet().
471: Input Parameters:
472: snes - the SNES context
473: its - iteration number
474: norm - 2-norm function value (may be estimated)
475: ctx - optional user-defined context for private data for the
476: monitor routine, as set by SNESMonitorSet()
478: Note:
479: See the manpage for PetscViewerDrawOpen() for useful runtime options,
480: such as -nox to deactivate all x-window output.
481: */
482: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
483: {
485: MonitorCtx *monP = (MonitorCtx*) ctx;
486: Vec x;
489: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %G\n",its,fnorm);
490: SNESGetSolution(snes,&x);
491: VecView(x,monP->viewer);
492: return(0);
493: }
495: /* ------------------------------------------------------------------- */
498: /*
499: PreCheck - Optional user-defined routine that checks the validity of
500: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
502: Input Parameters:
503: snes - the SNES context
504: ctx - optional user-defined context for private data for the
505: monitor routine, as set by SNESLineSearchSetPostCheck()
506: xcurrent - current solution
507: y - search direction and length
509: Output Parameters:
510: y - proposed step (search direction and length) (possibly changed)
511:
512: */
513: PetscErrorCode PreCheck(SNES snes,Vec xcurrent,Vec y,void *ctx,PetscTruth *changed_y)
514: {
516: *changed_y = PETSC_FALSE;
517: return(0);
518: }
520: /* ------------------------------------------------------------------- */
523: /*
524: PostCheck - Optional user-defined routine that checks the validity of
525: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
527: Input Parameters:
528: snes - the SNES context
529: ctx - optional user-defined context for private data for the
530: monitor routine, as set by SNESLineSearchSetPostCheck()
531: xcurrent - current solution
532: y - search direction and length
533: x - the new candidate iterate
535: Output Parameters:
536: y - proposed step (search direction and length) (possibly changed)
537: x - current iterate (possibly modified)
538:
539: */
540: PetscErrorCode PostCheck(SNES snes,Vec xcurrent,Vec y,Vec x,void *ctx,PetscTruth *changed_y,PetscTruth *changed_x)
541: {
543: PetscInt i,iter,xs,xm;
544: StepCheckCtx *check = (StepCheckCtx*) ctx;
545: ApplicationCtx *user = check->user;
546: PetscScalar *xa,*xa_last,tmp;
547: PetscReal rdiff;
548: DA da;
551: *changed_x = PETSC_FALSE;
552: *changed_y = PETSC_FALSE;
553: SNESGetIterationNumber(snes,&iter);
555: /* iteration 1 indicates we are working on the second iteration */
556: if (iter > 0) {
557: da = user->da;
558: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %G\n",iter,check->tolerance);
560: /* Access local array data */
561: DAVecGetArray(da,check->last_step,&xa_last);
562: DAVecGetArray(da,x,&xa);
563: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
565: /*
566: If we fail the user-defined check for validity of the candidate iterate,
567: then modify the iterate as we like. (Note that the particular modification
568: below is intended simply to demonstrate how to manipulate this data, not
569: as a meaningful or appropriate choice.)
570: */
571: for (i=xs; i<xs+xm; i++) {
572: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance; else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
573: if (rdiff > check->tolerance) {
574: tmp = xa[i];
575: xa[i] = .5*(xa[i] + xa_last[i]);
576: *changed_x = PETSC_TRUE;
577: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%G, x_last=%G, diff=%G, x_new=%G\n",
578: i,PetscAbsScalar(tmp),PetscAbsScalar(xa_last[i]),rdiff,PetscAbsScalar(xa[i]));
579: }
580: }
581: DAVecRestoreArray(da,check->last_step,&xa_last);
582: DAVecRestoreArray(da,x,&xa);
583: }
584: VecCopy(x,check->last_step);
585: return(0);
586: }