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