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;
35: PetscFunctionBeginUser;
36: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
37: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
38: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
39: PetscCall(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: PetscCall(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: PetscCall(VecCreate(comm, &x));
55: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
56: PetscCall(VecSetFromOptions(x));
57: PetscCall(VecDuplicate(x, &r));
58: PetscCall(VecDuplicate(x, &F));
59: PetscCall(VecDuplicate(x, &U));
61: /*
62: Set function evaluation routine and vector
63: */
64: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create matrix data structure; set Jacobian evaluation routine
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: PetscCall(MatCreate(comm, &J));
71: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
72: PetscCall(MatSetFromOptions(J));
73: PetscCall(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: PetscCall(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: PetscCall(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++) PetscCall(SNESConvergedReasonViewSet(snes, MySNESConvergedReasonView, &monP, 0));
100: PetscCall(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++) PetscCall(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: PetscCall(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: PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
120: v = xp * xp * xp;
121: PetscCall(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: PetscCall(FormInitialGuess(x));
135: PetscCall(SNESSolve(snes, NULL, x));
136: PetscCall(SNESGetIterationNumber(snes, &its));
138: /*
139: Free work space. All PETSc objects should be destroyed when they
140: are no longer needed.
141: */
142: PetscCall(VecDestroy(&x));
143: PetscCall(VecDestroy(&r));
144: PetscCall(VecDestroy(&U));
145: PetscCall(VecDestroy(&F));
146: PetscCall(MatDestroy(&J));
147: PetscCall(SNESDestroy(&snes));
148: PetscCall(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;
162: PetscFunctionBeginUser;
163: PetscCall(VecSet(x, pfive));
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
167: /*
168: FormFunction - Evaluates nonlinear function, F(x).
170: Input Parameters:
171: . snes - the SNES context
172: . x - input vector
173: . ctx - optional user-defined context, as set by SNESSetFunction()
175: Output Parameter:
176: . f - function vector
178: Note:
179: The user-defined context can contain any application-specific data
180: needed for the function evaluation (such as various parameters, work
181: vectors, and grid information). In this program the context is just
182: a vector containing the right-hand-side of the discretized PDE.
183: */
185: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
186: {
187: Vec g = (Vec)ctx;
188: const PetscScalar *xx, *gg;
189: PetscScalar *ff, d;
190: PetscInt i, n;
192: PetscFunctionBeginUser;
193: /*
194: Get pointers to vector data.
195: - For default PETSc vectors, VecGetArray() returns a pointer to
196: the data array. Otherwise, the routine is implementation dependent.
197: - You MUST call VecRestoreArray() when you no longer need access to
198: the array.
199: */
200: PetscCall(VecGetArrayRead(x, &xx));
201: PetscCall(VecGetArray(f, &ff));
202: PetscCall(VecGetArrayRead(g, &gg));
204: /*
205: Compute function
206: */
207: PetscCall(VecGetSize(x, &n));
208: d = (PetscReal)(n - 1);
209: d = d * d;
210: ff[0] = xx[0];
211: 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];
212: ff[n - 1] = xx[n - 1] - 1.0;
214: /*
215: Restore vectors
216: */
217: PetscCall(VecRestoreArrayRead(x, &xx));
218: PetscCall(VecRestoreArray(f, &ff));
219: PetscCall(VecRestoreArrayRead(g, &gg));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
222: /* ------------------------------------------------------------------- */
223: /*
224: FormJacobian - Evaluates Jacobian matrix.
226: Input Parameters:
227: . snes - the SNES context
228: . x - input vector
229: . dummy - optional user-defined context (not used here)
231: Output Parameters:
232: . jac - Jacobian matrix
233: . B - optionally different preconditioning matrix
235: */
237: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
238: {
239: const PetscScalar *xx;
240: PetscScalar A[3], d;
241: PetscInt i, n, j[3];
243: PetscFunctionBeginUser;
244: /*
245: Get pointer to vector data
246: */
247: PetscCall(VecGetArrayRead(x, &xx));
249: /*
250: Compute Jacobian entries and insert into matrix.
251: - Note that in this case we set all elements for a particular
252: row at once.
253: */
254: PetscCall(VecGetSize(x, &n));
255: d = (PetscReal)(n - 1);
256: d = d * d;
258: /*
259: Interior grid points
260: */
261: for (i = 1; i < n - 1; i++) {
262: j[0] = i - 1;
263: j[1] = i;
264: j[2] = i + 1;
265: A[0] = A[2] = d;
266: A[1] = -2.0 * d + 2.0 * xx[i];
267: PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
268: }
270: /*
271: Boundary points
272: */
273: i = 0;
274: A[0] = 1.0;
276: PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
278: i = n - 1;
279: A[0] = 1.0;
281: PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
283: /*
284: Restore vector
285: */
286: PetscCall(VecRestoreArrayRead(x, &xx));
288: /*
289: Assemble matrix
290: */
291: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
292: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
293: if (jac != B) {
294: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
295: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
296: }
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: PetscErrorCode MySNESConvergedReasonView(SNES snes, void *ctx)
301: {
302: ReasonViewCtx *monP = (ReasonViewCtx *)ctx;
303: PetscViewer viewer = monP->viewer;
304: SNESConvergedReason reason;
305: const char *strreason;
307: PetscFunctionBeginUser;
308: PetscCall(SNESGetConvergedReason(snes, &reason));
309: PetscCall(SNESGetConvergedReasonString(snes, &strreason));
310: PetscCall(PetscViewerASCIIPrintf(viewer, "Customized SNES converged reason view\n"));
311: PetscCall(PetscViewerASCIIAddTab(viewer, 1));
312: if (reason > 0) {
313: PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", strreason));
314: } else if (reason <= 0) {
315: PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", strreason));
316: }
317: PetscCall(PetscViewerASCIISubtractTab(viewer, 1));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: PetscErrorCode MyKSPConvergedReasonView(KSP ksp, void *ctx)
322: {
323: ReasonViewCtx *monP = (ReasonViewCtx *)ctx;
324: PetscViewer viewer = monP->viewer;
325: KSPConvergedReason reason;
326: const char *reasonstr;
328: PetscFunctionBeginUser;
329: PetscCall(KSPGetConvergedReason(ksp, &reason));
330: PetscCall(KSPGetConvergedReasonString(ksp, &reasonstr));
331: PetscCall(PetscViewerASCIIAddTab(viewer, 2));
332: PetscCall(PetscViewerASCIIPrintf(viewer, "Customized KSP converged reason view\n"));
333: PetscCall(PetscViewerASCIIAddTab(viewer, 1));
334: if (reason > 0) {
335: PetscCall(PetscViewerASCIIPrintf(viewer, "Converged due to %s\n", reasonstr));
336: } else if (reason <= 0) {
337: PetscCall(PetscViewerASCIIPrintf(viewer, "Did not converge due to %s\n", reasonstr));
338: }
339: PetscCall(PetscViewerASCIISubtractTab(viewer, 3));
340: PetscFunctionReturn(PETSC_SUCCESS);
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*/