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