Actual source code: ex6.c
2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
3: This example employs a user-defined reasonview routine.\n\n";
5: /*T
6: Concepts: SNES^basic uniprocessor example
7: Concepts: SNES^setting a user-defined reasonview routine
8: Processors: 1
9: T*/
11: #include <petscsnes.h>
13: /*
14: User-defined routines
15: */
16: PETSC_EXTERN PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
17: PETSC_EXTERN PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
18: PETSC_EXTERN PetscErrorCode FormInitialGuess(Vec);
19: PETSC_EXTERN PetscErrorCode MySNESConvergedReasonView(SNES,void*);
20: PETSC_EXTERN PetscErrorCode MyKSPConvergedReasonView(KSP,void*);
22: /*
23: User-defined context for monitoring
24: */
25: typedef struct {
26: PetscViewer viewer;
27: } ReasonViewCtx;
29: int main(int argc,char **argv)
30: {
31: SNES snes; /* SNES context */
32: KSP ksp; /* KSP context */
33: Vec x,r,F,U; /* vectors */
34: Mat J; /* Jacobian matrix */
35: ReasonViewCtx monP; /* monitoring context */
37: PetscInt its,n = 5,i;
38: PetscMPIInt size;
39: PetscScalar h,xp,v;
40: MPI_Comm comm;
42: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
43: MPI_Comm_size(PETSC_COMM_WORLD,&size);
44: if (size != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This is a uniprocessor example only!");
45: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
46: h = 1.0/(n-1);
47: comm = PETSC_COMM_WORLD;
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Create nonlinear solver context
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
52: SNESCreate(comm,&snes);
54: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
55: Create vector data structures; set function evaluation routine
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: /*
58: Note that we form 1 vector from scratch and then duplicate as needed.
59: */
60: VecCreate(comm,&x);
61: VecSetSizes(x,PETSC_DECIDE,n);
62: VecSetFromOptions(x);
63: VecDuplicate(x,&r);
64: VecDuplicate(x,&F);
65: VecDuplicate(x,&U);
67: /*
68: Set function evaluation routine and vector
69: */
70: SNESSetFunction(snes,r,FormFunction,(void*)F);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create matrix data structure; set Jacobian evaluation routine
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: MatCreate(comm,&J);
78: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
79: MatSetFromOptions(J);
80: MatSeqAIJSetPreallocation(J,3,NULL);
82: /*
83: Set Jacobian matrix data structure and default Jacobian evaluation
84: routine. User can override with:
85: -snes_fd : default finite differencing approximation of Jacobian
86: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
87: (unless user explicitly sets preconditioner)
88: -snes_mf_operator : form preconditioning matrix as set by the user,
89: but use matrix-free approx for Jacobian-vector
90: products within Newton-Krylov method
91: */
93: SNESSetJacobian(snes,J,J,FormJacobian,NULL);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Customize nonlinear solver; set runtime options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: /*
100: Set an optional user-defined reasonview routine
101: */
102: PetscViewerASCIIGetStdout(comm,&monP.viewer);
103: /* Just make sure we can not repeat addding the same function
104: * PETSc will be able to igore the repeated function
105: */
106: for (i=0; i<4; i++){
107: SNESConvergedReasonViewSet(snes,MySNESConvergedReasonView,&monP,0);
108: }
109: SNESGetKSP(snes,&ksp);
110: /* Just make sure we can not repeat addding the same function
111: * PETSc will be able to igore the repeated function
112: */
113: for (i=0; i<4; i++){
114: KSPConvergedReasonViewSet(ksp,MyKSPConvergedReasonView,&monP,0);
115: }
116: /*
117: Set SNES/KSP/KSP/PC runtime options, e.g.,
118: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
119: */
120: SNESSetFromOptions(snes);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Initialize application:
124: Store right-hand-side of PDE and exact solution
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: xp = 0.0;
128: for (i=0; i<n; i++) {
129: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
130: VecSetValues(F,1,&i,&v,INSERT_VALUES);
131: v = xp*xp*xp;
132: VecSetValues(U,1,&i,&v,INSERT_VALUES);
133: xp += h;
134: }
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Evaluate initial guess; then solve nonlinear system
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
139: /*
140: Note: The user should initialize the vector, x, with the initial guess
141: for the nonlinear solver prior to calling SNESSolve(). In particular,
142: to employ an initial guess of zero, the user should explicitly set
143: this vector to zero by calling VecSet().
144: */
145: FormInitialGuess(x);
146: SNESSolve(snes,NULL,x);
147: SNESGetIterationNumber(snes,&its);
149: /*
150: Free work space. All PETSc objects should be destroyed when they
151: are no longer needed.
152: */
153: VecDestroy(&x); VecDestroy(&r);
154: VecDestroy(&U); VecDestroy(&F);
155: MatDestroy(&J); SNESDestroy(&snes);
156: /*PetscViewerDestroy(&monP.viewer);*/
157: PetscFinalize();
158: return ierr;
159: }
160: /* ------------------------------------------------------------------- */
161: /*
162: FormInitialGuess - Computes initial guess.
164: Input/Output Parameter:
165: . x - the solution vector
166: */
167: PetscErrorCode FormInitialGuess(Vec x)
168: {
170: PetscScalar pfive = .50;
171: VecSet(x,pfive);
172: return 0;
173: }
174: /* ------------------------------------------------------------------- */
175: /*
176: FormFunction - Evaluates nonlinear function, F(x).
178: Input Parameters:
179: . snes - the SNES context
180: . x - input vector
181: . ctx - optional user-defined context, as set by SNESSetFunction()
183: Output Parameter:
184: . f - function vector
186: Note:
187: The user-defined context can contain any application-specific data
188: needed for the function evaluation (such as various parameters, work
189: vectors, and grid information). In this program the context is just
190: a vector containing the right-hand-side of the discretized PDE.
191: */
193: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
194: {
195: Vec g = (Vec)ctx;
196: const PetscScalar *xx,*gg;
197: PetscScalar *ff,d;
198: PetscErrorCode ierr;
199: PetscInt i,n;
201: /*
202: Get pointers to vector data.
203: - For default PETSc vectors, VecGetArray() returns a pointer to
204: the data array. Otherwise, the routine is implementation dependent.
205: - You MUST call VecRestoreArray() when you no longer need access to
206: the array.
207: */
208: VecGetArrayRead(x,&xx);
209: VecGetArray(f,&ff);
210: VecGetArrayRead(g,&gg);
212: /*
213: Compute function
214: */
215: VecGetSize(x,&n);
216: d = (PetscReal)(n - 1); d = d*d;
217: ff[0] = xx[0];
218: for (i=1; i<n-1; i++) ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
219: ff[n-1] = xx[n-1] - 1.0;
221: /*
222: Restore vectors
223: */
224: VecRestoreArrayRead(x,&xx);
225: VecRestoreArray(f,&ff);
226: VecRestoreArrayRead(g,&gg);
227: return 0;
228: }
229: /* ------------------------------------------------------------------- */
230: /*
231: FormJacobian - Evaluates Jacobian matrix.
233: Input Parameters:
234: . snes - the SNES context
235: . x - input vector
236: . dummy - optional user-defined context (not used here)
238: Output Parameters:
239: . jac - Jacobian matrix
240: . B - optionally different preconditioning matrix
242: */
244: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
245: {
246: const PetscScalar *xx;
247: PetscScalar A[3],d;
248: PetscErrorCode ierr;
249: PetscInt i,n,j[3];
251: /*
252: Get pointer to vector data
253: */
254: VecGetArrayRead(x,&xx);
256: /*
257: Compute Jacobian entries and insert into matrix.
258: - Note that in this case we set all elements for a particular
259: row at once.
260: */
261: VecGetSize(x,&n);
262: d = (PetscReal)(n - 1); d = d*d;
264: /*
265: Interior grid points
266: */
267: for (i=1; i<n-1; i++) {
268: j[0] = i - 1; j[1] = i; j[2] = i + 1;
269: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
270: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
271: }
273: /*
274: Boundary points
275: */
276: i = 0; A[0] = 1.0;
278: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
280: i = n-1; A[0] = 1.0;
282: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
284: /*
285: Restore vector
286: */
287: VecRestoreArrayRead(x,&xx);
289: /*
290: Assemble matrix
291: */
292: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
293: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
294: if (jac != B) {
295: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
296: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
297: }
298: return 0;
299: }
301: PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx)
302: {
303: PetscErrorCode ierr;
304: ReasonViewCtx *monP = (ReasonViewCtx*) ctx;
305: PetscViewer viewer = monP->viewer;
306: SNESConvergedReason reason;
307: const char *strreason;
309: SNESGetConvergedReason(snes,&reason);
310: SNESGetConvergedReasonString(snes,&strreason);
311: PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");
312: PetscViewerASCIIAddTab(viewer,1);
313: if (reason > 0) {
314: PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);
315: } else if (reason <= 0) {
316: PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);
317: }
318: PetscViewerASCIISubtractTab(viewer,1);
319: return 0;
320: }
322: PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx)
323: {
324: PetscErrorCode ierr;
325: ReasonViewCtx *monP = (ReasonViewCtx*) ctx;
326: PetscViewer viewer = monP->viewer;
327: KSPConvergedReason reason;
328: const char *reasonstr;
330: KSPGetConvergedReason(ksp,&reason);
331: KSPGetConvergedReasonString(ksp,&reasonstr);
332: PetscViewerASCIIAddTab(viewer,2);
333: PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n");
334: PetscViewerASCIIAddTab(viewer,1);
335: if (reason > 0) {
336: PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr);
337: } else if (reason <= 0) {
338: PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr);
339: }
340: PetscViewerASCIISubtractTab(viewer,3);
341: return 0;
342: }
345: /*TEST
347: test:
348: suffix: 1
349: nsize: 1
351: test:
352: suffix: 2
353: nsize: 1
354: args: -ksp_converged_reason_view_cancel
356: test:
357: suffix: 3
358: nsize: 1
359: args: -ksp_converged_reason_view_cancel -ksp_converged_reason
361: test:
362: suffix: 4
363: nsize: 1
364: args: -snes_converged_reason_view_cancel
366: test:
367: suffix: 5
368: nsize: 1
369: args: -snes_converged_reason_view_cancel -snes_converged_reason
371: TEST*/