Actual source code: ex3.c
1: /*$Id: ex3.c,v 1.87 2001/08/07 21:31:17 bsmith Exp $*/
3: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.n
4: This example employs a user-defined monitoring routine and optionally a user-definedn
5: routine to check candidate iterates produced by line search routines. This code alson
6: demonstrates use of the macro __FUNCT__ to define routine names for use in error handling.n
7: The command line options include:n
8: -check_iterates : activate checking of iteratesn
9: -check_tol <tol>: set tolerance for iterate checkingnn";
11: /*T
12: Concepts: SNES^basic parallel example
13: Concepts: SNES^setting a user-defined monitoring routine
14: Concepts: error handling^using the macro __FUNCT__ to define routine names;
15: Processors: n
16: T*/
18: /*
19: Include "petscdraw.h" so that we can use distributed arrays (DAs).
20: Include "petscdraw.h" so that we can use PETSc drawing routines.
21: Include "petscsnes.h" so that we can use SNES solvers. Note that this
22: file automatically includes:
23: petsc.h - base PETSc routines petscvec.h - vectors
24: petscsys.h - system routines petscmat.h - matrices
25: petscis.h - index sets petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers petscpc.h - preconditioners
27: petscsles.h - linear solvers
28: */
30: #include petscda.h
31: #include petscsnes.h
33: /*
34: User-defined routines. Note that immediately before each routine below,
35: we define the macro __FUNCT__ to be a string containing the routine name.
36: If defined, this macro is used in the PETSc error handlers to provide a
37: complete traceback of routine names. All PETSc library routines use this
38: macro, and users can optionally employ it as well in their application
39: codes. Note that users can get a traceback of PETSc errors regardless of
40: whether they define __FUNCT__ in application codes; this macro merely
41: provides the added traceback detail of the application routine names.
42: */
43: int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
44: int FormFunction(SNES,Vec,Vec,void*);
45: int FormInitialGuess(Vec);
46: int Monitor(SNES,int,PetscReal,void *);
47: int StepCheck(SNES,void *,Vec,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: int rank; /* rank of processor */
56: int 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: } StepCheckCtx;
76: #undef __FUNCT__
78: int main(int argc,char **argv)
79: {
80: SNES snes; /* SNES context */
81: Mat J; /* Jacobian matrix */
82: ApplicationCtx ctx; /* user-defined context */
83: Vec x,r,U,F; /* vectors */
84: MonitorCtx monP; /* monitoring context */
85: StepCheckCtx checkP; /* step-checking context */
86: PetscTruth step_check; /* flag indicating whether we're checking
87: candidate iterates */
88: PetscScalar xp,*FF,*UU,none = -1.0;
89: int ierr,its,N = 5,i,maxit,maxf,xs,xm;
90: PetscReal atol,rtol,stol,norm;
92: PetscInitialize(&argc,&argv,(char *)0,help);
93: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
94: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
95: PetscOptionsGetInt(PETSC_NULL,"-n",&N,PETSC_NULL);
96: ctx.h = 1.0/(N-1);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Create nonlinear solver context
100: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: SNESCreate(PETSC_COMM_WORLD,&snes);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Create vector data structures; set function evaluation routine
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: /*
109: Create distributed array (DA) to manage parallel grid and vectors
110: */
111: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,N,1,1,PETSC_NULL,&ctx.da);
113: /*
114: Extract global and local vectors from DA; then duplicate for remaining
115: vectors that are the same types
116: */
117: DACreateGlobalVector(ctx.da,&x);
118: VecDuplicate(x,&r);
119: VecDuplicate(x,&F); ctx.F = F;
120: VecDuplicate(x,&U);
122: /*
123: Set function evaluation routine and vector. Whenever the nonlinear
124: solver needs to compute the nonlinear function, it will call this
125: routine.
126: - Note that the final routine argument is the user-defined
127: context that provides application-specific data for the
128: function evaluation routine.
129: */
130: SNESSetFunction(snes,r,FormFunction,&ctx);
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Create matrix data structure; set Jacobian evaluation routine
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&J);
137: MatSetFromOptions(J);
139: /*
140: Set Jacobian matrix data structure and default Jacobian evaluation
141: routine. Whenever the nonlinear solver needs to compute the
142: Jacobian matrix, it will call this routine.
143: - Note that the final routine argument is the user-defined
144: context that provides application-specific data for the
145: Jacobian evaluation routine.
146: */
147: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
149: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150: Customize nonlinear solver; set runtime options
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: /*
154: Set an optional user-defined monitoring routine
155: */
156: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
157: SNESSetMonitor(snes,Monitor,&monP,0);
159: /*
160: Set names for some vectors to facilitate monitoring (optional)
161: */
162: PetscObjectSetName((PetscObject)x,"Approximate Solution");
163: PetscObjectSetName((PetscObject)U,"Exact Solution");
165: /*
166: Set SNES/SLES/KSP/PC runtime options, e.g.,
167: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
168: */
169: SNESSetFromOptions(snes);
171: /*
172: Set an optional user-defined routine to check the validity of candidate
173: iterates that are determined by line search methods
174: */
175: PetscOptionsHasName(PETSC_NULL,"-check_iterates",&step_check);
176: if (step_check) {
177: PetscPrintf(PETSC_COMM_WORLD,"Activating step checking routinen");
178: SNESSetLineSearchCheck(snes,StepCheck,&checkP);
179: VecDuplicate(x,&(checkP.last_step));
180: checkP.tolerance = 1.0;
181: PetscOptionsGetReal(PETSC_NULL,"-check_tol",&checkP.tolerance,PETSC_NULL);
182: }
185: /*
186: Print parameters used for convergence testing (optional) ... just
187: to demonstrate this routine; this information is also printed with
188: the option -snes_view
189: */
190: SNESGetTolerances(snes,&atol,&rtol,&stol,&maxit,&maxf);
191: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%d, maxf=%dn",atol,rtol,stol,maxit,maxf);
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Initialize application:
195: Store right-hand-side of PDE and exact solution
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: /*
199: Get local grid boundaries (for 1-dimensional DA):
200: xs, xm - starting grid index, width of local grid (no ghost points)
201: */
202: DAGetCorners(ctx.da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
204: /*
205: Get pointers to vector data
206: */
207: DAVecGetArray(ctx.da,F,(void**)&FF);
208: DAVecGetArray(ctx.da,U,(void**)&UU);
210: /*
211: Compute local vector entries
212: */
213: xp = ctx.h*xs;
214: for (i=xs; i<xs+xm; i++) {
215: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
216: UU[i] = xp*xp*xp;
217: xp += ctx.h;
218: }
220: /*
221: Restore vectors
222: */
223: DAVecRestoreArray(ctx.da,F,(void**)&FF);
224: DAVecRestoreArray(ctx.da,U,(void**)&UU);
226: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227: Evaluate initial guess; then solve nonlinear system
228: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230: /*
231: Note: The user should initialize the vector, x, with the initial guess
232: for the nonlinear solver prior to calling SNESSolve(). In particular,
233: to employ an initial guess of zero, the user should explicitly set
234: this vector to zero by calling VecSet().
235: */
236: FormInitialGuess(x);
237: SNESSolve(snes,x,&its);
238: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %dnn",its);
240: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: Check solution and clean up
242: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
244: /*
245: Check the error
246: */
247: VecAXPY(&none,U,x);
248: ierr = VecNorm(x,NORM_2,&norm);
249: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A, Iterations %dn",norm,its);
251: /*
252: Free work space. All PETSc objects should be destroyed when they
253: are no longer needed.
254: */
255: PetscViewerDestroy(monP.viewer);
256: if (step_check) {VecDestroy(checkP.last_step);}
257: VecDestroy(x);
258: VecDestroy(r);
259: VecDestroy(U);
260: VecDestroy(F);
261: MatDestroy(J);
262: SNESDestroy(snes);
263: DADestroy(ctx.da);
264: PetscFinalize();
266: return(0);
267: }
268: /* ------------------------------------------------------------------- */
269: #undef __FUNCT__
271: /*
272: FormInitialGuess - Computes initial guess.
274: Input/Output Parameter:
275: . x - the solution vector
276: */
277: int FormInitialGuess(Vec x)
278: {
279: int ierr;
280: PetscScalar pfive = .50;
283: VecSet(&pfive,x);
284: return(0);
285: }
286: /* ------------------------------------------------------------------- */
287: #undef __FUNCT__
289: /*
290: FormFunction - Evaluates nonlinear function, F(x).
292: Input Parameters:
293: . snes - the SNES context
294: . x - input vector
295: . ctx - optional user-defined context, as set by SNESSetFunction()
297: Output Parameter:
298: . f - function vector
300: Note:
301: The user-defined context can contain any application-specific
302: data needed for the function evaluation.
303: */
304: int FormFunction(SNES snes,Vec x,Vec f,void *ctx)
305: {
306: ApplicationCtx *user = (ApplicationCtx*) ctx;
307: DA da = user->da;
308: PetscScalar *xx,*ff,*FF,d;
309: int i,ierr,M,xs,xm;
310: Vec xlocal;
313: DAGetLocalVector(da,&xlocal);
314: /*
315: Scatter ghost points to local vector, using the 2-step process
316: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
317: By placing code between these two statements, computations can
318: be done while messages are in transition.
319: */
320: DAGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
321: DAGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
323: /*
324: Get pointers to vector data.
325: - The vector xlocal includes ghost point; the vectors x and f do
326: NOT include ghost points.
327: - Using DAVecGetArray() allows accessing the values using global ordering
328: */
329: DAVecGetArray(da,xlocal,(void**)&xx);
330: DAVecGetArray(da,f,(void**)&ff);
331: DAVecGetArray(da,user->F,(void**)&FF);
333: /*
334: Get local grid boundaries (for 1-dimensional DA):
335: xs, xm - starting grid index, width of local grid (no ghost points)
336: */
337: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
338: DAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
339: PETSC_NULL,PETSC_NULL,PETSC_NULL);
341: /*
342: Set function values for boundary points; define local interior grid point range:
343: xsi - starting interior grid index
344: xei - ending interior grid index
345: */
346: if (xs == 0) { /* left boundary */
347: ff[0] = xx[0];
348: xs++;xm--;
349: }
350: if (xs+xm == M) { /* right boundary */
351: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
352: xm--;
353: }
355: /*
356: Compute function over locally owned part of the grid (interior points only)
357: */
358: d = 1.0/(user->h*user->h);
359: for (i=xs; i<xs+xm; i++) {
360: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
361: }
363: /*
364: Restore vectors
365: */
366: DAVecRestoreArray(da,xlocal,(void**)&xx);
367: DAVecRestoreArray(da,f,(void**)&ff);
368: DAVecRestoreArray(da,user->F,(void**)&FF);
369: DARestoreLocalVector(da,&xlocal);
370: return(0);
371: }
372: /* ------------------------------------------------------------------- */
373: #undef __FUNCT__
375: /*
376: FormJacobian - Evaluates Jacobian matrix.
378: Input Parameters:
379: . snes - the SNES context
380: . x - input vector
381: . dummy - optional user-defined context (not used here)
383: Output Parameters:
384: . jac - Jacobian matrix
385: . B - optionally different preconditioning matrix
386: . flag - flag indicating matrix structure
387: */
388: int FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *ctx)
389: {
390: ApplicationCtx *user = (ApplicationCtx*) ctx;
391: PetscScalar *xx,d,A[3];
392: int i,j[3],ierr,M,xs,xm;
393: DA da = user->da;
396: /*
397: Get pointer to vector data
398: */
399: DAVecGetArray(da,x,(void**)&xx);
400: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
402: /*
403: Get range of locally owned matrix
404: */
405: DAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
406: PETSC_NULL,PETSC_NULL,PETSC_NULL);
408: /*
409: Determine starting and ending local indices for interior grid points.
410: Set Jacobian entries for boundary points.
411: */
413: if (xs == 0) { /* left boundary */
414: i = 0; A[0] = 1.0;
415: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
416: xs++;xm--;
417: }
418: if (xs+xm == M) { /* right boundary */
419: i = M-1; A[0] = 1.0;
420: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
421: xm--;
422: }
424: /*
425: Interior grid points
426: - Note that in this case we set all elements for a particular
427: row at once.
428: */
429: d = 1.0/(user->h*user->h);
430: for (i=xs; i<xs+xm; i++) {
431: j[0] = i - 1; j[1] = i; j[2] = i + 1;
432: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
433: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
434: }
436: /*
437: Assemble matrix, using the 2-step process:
438: MatAssemblyBegin(), MatAssemblyEnd().
439: By placing code between these two statements, computations can be
440: done while messages are in transition.
442: Also, restore vector.
443: */
445: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
446: DAVecRestoreArray(da,x,(void**)&xx);
447: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
449: *flag = SAME_NONZERO_PATTERN;
450: return(0);
451: }
452: /* ------------------------------------------------------------------- */
453: #undef __FUNCT__
455: /*
456: Monitor - Optional user-defined monitoring routine that views the
457: current iterate with an x-window plot. Set by SNESSetMonitor().
459: Input Parameters:
460: snes - the SNES context
461: its - iteration number
462: norm - 2-norm function value (may be estimated)
463: ctx - optional user-defined context for private data for the
464: monitor routine, as set by SNESSetMonitor()
466: Note:
467: See the manpage for PetscViewerDrawOpen() for useful runtime options,
468: such as -nox to deactivate all x-window output.
469: */
470: int Monitor(SNES snes,int its,PetscReal fnorm,void *ctx)
471: {
472: int ierr;
473: MonitorCtx *monP = (MonitorCtx*) ctx;
474: Vec x;
477: PetscPrintf(PETSC_COMM_WORLD,"iter = %d,SNES Function norm %gn",its,fnorm);
478: SNESGetSolution(snes,&x);
479: VecView(x,monP->viewer);
480: return(0);
481: }
482: /* ------------------------------------------------------------------- */
483: #undef __FUNCT__
485: /*
486: StepCheck - Optional user-defined routine that checks the validity of
487: candidate steps of a line search method. Set by SNESSetLineSearchCheck().
489: Input Parameters:
490: snes - the SNES context
491: ctx - optional user-defined context for private data for the
492: monitor routine, as set by SNESSetLineSearchCheck()
493: x - the new candidate iterate
495: Output Parameters:
496: x - current iterate (possibly modified)
497: flg - flag indicating whether x has been modified (either
498: PETSC_TRUE of PETSC_FALSE)
499: */
500: int StepCheck(SNES snes,void *ctx,Vec x,PetscTruth *flg)
501: {
502: int ierr,i,iter,xs,xm;
503: ApplicationCtx *user;
504: StepCheckCtx *check = (StepCheckCtx*) ctx;
505: PetscScalar *xa,*xa_last,tmp;
506: PetscReal rdiff;
507: DA da;
510: *flg = PETSC_FALSE;
511: SNESGetIterationNumber(snes,&iter);
513: if (iter > 1) {
514: SNESGetApplicationContext(snes,(void**)&user);
515: da = user->da;
516: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %d with tolerance %gn",
517: iter,check->tolerance);
519: /* Access local array data */
520: DAVecGetArray(da,check->last_step,(void**)&xa_last);
521: DAVecGetArray(da,x,(void**)&xa);
522: DAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
524: /*
525: If we fail the user-defined check for validity of the candidate iterate,
526: then modify the iterate as we like. (Note that the particular modification
527: below is intended simply to demonstrate how to manipulate this data, not
528: as a meaningful or appropriate choice.)
529: */
530: for (i=xs; i<xs+xm; i++) {
531: rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
532: if (rdiff > check->tolerance) {
533: tmp = xa[i];
534: xa[i] = (xa[i] + xa_last[i])/2.0;
535: *flg = PETSC_TRUE;
536: PetscPrintf(PETSC_COMM_WORLD," Altering entry %d: x=%g, x_last=%g, diff=%g, x_new=%gn",
537: i,PetscAbsScalar(tmp),PetscAbsScalar(xa_last[i]),rdiff,PetscAbsScalar(xa[i]));
538: }
539: }
540: DAVecRestoreArray(da,check->last_step,(void**)&xa_last);
541: DAVecRestoreArray(da,x,(void**)&xa);
542: }
543: VecCopy(x,check->last_step);
545: return(0);
546: }