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);
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create matrix data structure; set Jacobian evaluation routine
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: MatCreate(comm,&J);
77: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);
78: MatSetFromOptions(J);
79: MatSeqAIJSetPreallocation(J,3,NULL);
81: /*
82: Set Jacobian matrix data structure and default Jacobian evaluation
83: routine. User can override with:
84: -snes_fd : default finite differencing approximation of Jacobian
85: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
86: (unless user explicitly sets preconditioner)
87: -snes_mf_operator : form preconditioning matrix as set by the user,
88: but use matrix-free approx for Jacobian-vector
89: products within Newton-Krylov method
90: */
92: SNESSetJacobian(snes,J,J,FormJacobian,NULL);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Customize nonlinear solver; set runtime options
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: /*
99: Set an optional user-defined reasonview routine
100: */
101: PetscViewerASCIIGetStdout(comm,&monP.viewer);
102: /* Just make sure we can not repeat addding the same function
103: * PETSc will be able to igore the repeated function
104: */
105: for (i=0; i<4; i++) {
106: SNESConvergedReasonViewSet(snes,MySNESConvergedReasonView,&monP,0);
107: }
108: SNESGetKSP(snes,&ksp);
109: /* Just make sure we can not repeat addding the same function
110: * PETSc will be able to igore the repeated function
111: */
112: for (i=0; i<4; i++) {
113: KSPConvergedReasonViewSet(ksp,MyKSPConvergedReasonView,&monP,0);
114: }
115: /*
116: Set SNES/KSP/KSP/PC runtime options, e.g.,
117: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
118: */
119: SNESSetFromOptions(snes);
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Initialize application:
123: Store right-hand-side of PDE and exact solution
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: xp = 0.0;
127: for (i=0; i<n; i++) {
128: v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
129: VecSetValues(F,1,&i,&v,INSERT_VALUES);
130: v = xp*xp*xp;
131: VecSetValues(U,1,&i,&v,INSERT_VALUES);
132: xp += h;
133: }
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Evaluate initial guess; then solve nonlinear system
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: /*
139: Note: The user should initialize the vector, x, with the initial guess
140: for the nonlinear solver prior to calling SNESSolve(). In particular,
141: to employ an initial guess of zero, the user should explicitly set
142: this vector to zero by calling VecSet().
143: */
144: FormInitialGuess(x);
145: SNESSolve(snes,NULL,x);
146: SNESGetIterationNumber(snes,&its);
148: /*
149: Free work space. All PETSc objects should be destroyed when they
150: are no longer needed.
151: */
152: VecDestroy(&x); VecDestroy(&r);
153: VecDestroy(&U); VecDestroy(&F);
154: MatDestroy(&J); SNESDestroy(&snes);
155: /*PetscViewerDestroy(&monP.viewer);*/
156: PetscFinalize();
157: return ierr;
158: }
159: /* ------------------------------------------------------------------- */
160: /*
161: FormInitialGuess - Computes initial guess.
163: Input/Output Parameter:
164: . x - the solution vector
165: */
166: PetscErrorCode FormInitialGuess(Vec x)
167: {
169: PetscScalar pfive = .50;
170: VecSet(x,pfive);
171: return 0;
172: }
173: /* ------------------------------------------------------------------- */
174: /*
175: FormFunction - Evaluates nonlinear function, F(x).
177: Input Parameters:
178: . snes - the SNES context
179: . x - input vector
180: . ctx - optional user-defined context, as set by SNESSetFunction()
182: Output Parameter:
183: . f - function vector
185: Note:
186: The user-defined context can contain any application-specific data
187: needed for the function evaluation (such as various parameters, work
188: vectors, and grid information). In this program the context is just
189: a vector containing the right-hand-side of the discretized PDE.
190: */
192: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
193: {
194: Vec g = (Vec)ctx;
195: const PetscScalar *xx,*gg;
196: PetscScalar *ff,d;
197: PetscErrorCode ierr;
198: PetscInt i,n;
200: /*
201: Get pointers to vector data.
202: - For default PETSc vectors, VecGetArray() returns a pointer to
203: the data array. Otherwise, the routine is implementation dependent.
204: - You MUST call VecRestoreArray() when you no longer need access to
205: the array.
206: */
207: VecGetArrayRead(x,&xx);
208: VecGetArray(f,&ff);
209: VecGetArrayRead(g,&gg);
211: /*
212: Compute function
213: */
214: VecGetSize(x,&n);
215: d = (PetscReal)(n - 1); d = d*d;
216: ff[0] = xx[0];
217: 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];
218: ff[n-1] = xx[n-1] - 1.0;
220: /*
221: Restore vectors
222: */
223: VecRestoreArrayRead(x,&xx);
224: VecRestoreArray(f,&ff);
225: VecRestoreArrayRead(g,&gg);
226: return 0;
227: }
228: /* ------------------------------------------------------------------- */
229: /*
230: FormJacobian - Evaluates Jacobian matrix.
232: Input Parameters:
233: . snes - the SNES context
234: . x - input vector
235: . dummy - optional user-defined context (not used here)
237: Output Parameters:
238: . jac - Jacobian matrix
239: . B - optionally different preconditioning matrix
241: */
243: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
244: {
245: const PetscScalar *xx;
246: PetscScalar A[3],d;
247: PetscErrorCode ierr;
248: PetscInt i,n,j[3];
250: /*
251: Get pointer to vector data
252: */
253: VecGetArrayRead(x,&xx);
255: /*
256: Compute Jacobian entries and insert into matrix.
257: - Note that in this case we set all elements for a particular
258: row at once.
259: */
260: VecGetSize(x,&n);
261: d = (PetscReal)(n - 1); d = d*d;
263: /*
264: Interior grid points
265: */
266: for (i=1; i<n-1; i++) {
267: j[0] = i - 1; j[1] = i; j[2] = i + 1;
268: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
269: MatSetValues(B,1,&i,3,j,A,INSERT_VALUES);
270: }
272: /*
273: Boundary points
274: */
275: i = 0; A[0] = 1.0;
277: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
279: i = n-1; A[0] = 1.0;
281: MatSetValues(B,1,&i,1,&i,A,INSERT_VALUES);
283: /*
284: Restore vector
285: */
286: VecRestoreArrayRead(x,&xx);
288: /*
289: Assemble matrix
290: */
291: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
292: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
293: if (jac != B) {
294: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
295: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
296: }
297: return 0;
298: }
300: PetscErrorCode MySNESConvergedReasonView(SNES snes,void *ctx)
301: {
302: PetscErrorCode ierr;
303: ReasonViewCtx *monP = (ReasonViewCtx*) ctx;
304: PetscViewer viewer = monP->viewer;
305: SNESConvergedReason reason;
306: const char *strreason;
308: SNESGetConvergedReason(snes,&reason);
309: SNESGetConvergedReasonString(snes,&strreason);
310: PetscViewerASCIIPrintf(viewer,"Customized SNES converged reason view\n");
311: PetscViewerASCIIAddTab(viewer,1);
312: if (reason > 0) {
313: PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",strreason);
314: } else if (reason <= 0) {
315: PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",strreason);
316: }
317: PetscViewerASCIISubtractTab(viewer,1);
318: return 0;
319: }
321: PetscErrorCode MyKSPConvergedReasonView(KSP ksp,void *ctx)
322: {
323: PetscErrorCode ierr;
324: ReasonViewCtx *monP = (ReasonViewCtx*) ctx;
325: PetscViewer viewer = monP->viewer;
326: KSPConvergedReason reason;
327: const char *reasonstr;
329: KSPGetConvergedReason(ksp,&reason);
330: KSPGetConvergedReasonString(ksp,&reasonstr);
331: PetscViewerASCIIAddTab(viewer,2);
332: PetscViewerASCIIPrintf(viewer,"Customized KSP converged reason view\n");
333: PetscViewerASCIIAddTab(viewer,1);
334: if (reason > 0) {
335: PetscViewerASCIIPrintf(viewer,"Converged due to %s\n",reasonstr);
336: } else if (reason <= 0) {
337: PetscViewerASCIIPrintf(viewer,"Did not converge due to %s\n",reasonstr);
338: }
339: PetscViewerASCIISubtractTab(viewer,3);
340: return 0;
341: }
343: /*TEST
345: test:
346: suffix: 1
347: nsize: 1
348: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"
350: test:
351: suffix: 2
352: nsize: 1
353: args: -ksp_converged_reason_view_cancel
354: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"
356: test:
357: suffix: 3
358: nsize: 1
359: args: -ksp_converged_reason_view_cancel -ksp_converged_reason
360: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"
362: test:
363: suffix: 4
364: nsize: 1
365: args: -snes_converged_reason_view_cancel
366: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"
368: test:
369: suffix: 5
370: nsize: 1
371: args: -snes_converged_reason_view_cancel -snes_converged_reason
372: filter: sed -e "s/CONVERGED_ATOL/CONVERGED_RTOL/g"
374: TEST*/