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";
3
4 #include <petscsnes.h>
5
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 *);
14
15 /*
16 User-defined context for monitoring
17 */
18 typedef struct {
19 PetscViewer viewer;
20 } ReasonViewCtx;
21
main(int argc,char ** argv)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;
33
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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44
45 PetscCall(SNESCreate(comm, &snes));
46
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));
59
60 /*
61 Set function evaluation routine and vector
62 */
63 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
64
65 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66 Create matrix data structure; set Jacobian evaluation routine
67 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68
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));
73
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 matrix used to construct the preconditioner as set by the user,
81 but use matrix-free approx for Jacobian-vector
82 products within Newton-Krylov method
83 */
84
85 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
86
87 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88 Customize nonlinear solver; set runtime options
89 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90
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));
109
110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111 Initialize application:
112 Store right-hand side of PDE and exact solution
113 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114
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 }
123
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));
136
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 }
150
151 /*
152 FormInitialGuess - Computes initial guess.
153
154 Input/Output Parameter:
155 . x - the solution vector
156 */
FormInitialGuess(Vec x)157 PetscErrorCode FormInitialGuess(Vec x)
158 {
159 PetscScalar pfive = .50;
160
161 PetscFunctionBeginUser;
162 PetscCall(VecSet(x, pfive));
163 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165
166 /*
167 FormFunction - Evaluates nonlinear function, F(x).
168
169 Input Parameters:
170 . snes - the SNES context
171 . x - input vector
172 . ctx - optional user-defined context, as set by SNESSetFunction()
173
174 Output Parameter:
175 . f - function vector
176
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 */
183
FormFunction(SNES snes,Vec x,Vec f,PetscCtx ctx)184 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
185 {
186 Vec g = (Vec)ctx;
187 const PetscScalar *xx, *gg;
188 PetscScalar *ff, d;
189 PetscInt i, n;
190
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));
202
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;
212
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.
224
225 Input Parameters:
226 . snes - the SNES context
227 . x - input vector
228 . dummy - optional user-defined context (not used here)
229
230 Output Parameters:
231 . jac - Jacobian matrix
232 . B - optionally different matrix used to construct the preconditioner
233
234 */
235
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)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];
241
242 PetscFunctionBeginUser;
243 /*
244 Get pointer to vector data
245 */
246 PetscCall(VecGetArrayRead(x, &xx));
247
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;
256
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 }
268
269 /*
270 Boundary points
271 */
272 i = 0;
273 A[0] = 1.0;
274
275 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
276
277 i = n - 1;
278 A[0] = 1.0;
279
280 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
281
282 /*
283 Restore vector
284 */
285 PetscCall(VecRestoreArrayRead(x, &xx));
286
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 }
298
MySNESConvergedReasonView(SNES snes,PetscCtx ctx)299 PetscErrorCode MySNESConvergedReasonView(SNES snes, PetscCtx ctx)
300 {
301 ReasonViewCtx *monP = (ReasonViewCtx *)ctx;
302 PetscViewer viewer = monP->viewer;
303 SNESConvergedReason reason;
304 const char *strreason;
305
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 }
319
MyKSPConvergedReasonView(KSP ksp,PetscCtx ctx)320 PetscErrorCode MyKSPConvergedReasonView(KSP ksp, PetscCtx ctx)
321 {
322 ReasonViewCtx *monP = (ReasonViewCtx *)ctx;
323 PetscViewer viewer = monP->viewer;
324 KSPConvergedReason reason;
325 const char *reasonstr;
326
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 }
341
342 /*TEST
343
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"
348
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"
354
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"
360
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"
366
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"
372
373 TEST*/
374