1 static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
2 This example employs a user-defined monitoring routine.\n\n";
3
4 /*
5 Include "petscdraw.h" so that we can use PETSc drawing routines.
6 Include "petscsnes.h" so that we can use SNES solvers. Note that this
7 file automatically includes:
8 petscsys.h - base PETSc routines petscvec.h - vectors
9 petscmat.h - matrices
10 petscis.h - index sets petscksp.h - Krylov subspace methods
11 petscviewer.h - viewers petscpc.h - preconditioners
12 petscksp.h - linear solvers
13 */
14
15 #include <petscsnes.h>
16
17 /*
18 User-defined routines
19 */
20 extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
21 extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
22 extern PetscErrorCode FormInitialGuess(Vec);
23 extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
24
25 /*
26 User-defined context for monitoring
27 */
28 typedef struct {
29 PetscViewer viewer;
30 } MonitorCtx;
31
main(int argc,char ** argv)32 int main(int argc, char **argv)
33 {
34 SNES snes; /* SNES context */
35 Vec x, r, F, U; /* vectors */
36 Mat J; /* Jacobian matrix */
37 MonitorCtx monP; /* monitoring context */
38 PetscInt its, n = 5, i, maxit, maxf;
39 PetscMPIInt size;
40 PetscScalar h, xp, v, none = -1.0;
41 PetscReal abstol, rtol, stol, norm;
42
43 PetscFunctionBeginUser;
44 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
45 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
46 PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
47 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
48 h = 1.0 / (n - 1);
49
50 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51 Create nonlinear solver context
52 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53
54 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
55
56 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57 Create vector data structures; set function evaluation routine
58 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59 /*
60 Note that we form 1 vector from scratch and then duplicate as needed.
61 */
62 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
63 PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
64 PetscCall(VecSetFromOptions(x));
65 PetscCall(VecDuplicate(x, &r));
66 PetscCall(VecDuplicate(x, &F));
67 PetscCall(VecDuplicate(x, &U));
68
69 /*
70 Set function evaluation routine and vector
71 */
72 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
73
74 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75 Create matrix data structure; set Jacobian evaluation routine
76 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77
78 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
79 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n));
80 PetscCall(MatSetFromOptions(J));
81 PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
82
83 /*
84 Set Jacobian matrix data structure and default Jacobian evaluation
85 routine. User can override with:
86 -snes_fd : default finite differencing approximation of Jacobian
87 -snes_mf : matrix-free Newton-Krylov method with no preconditioning
88 (unless user explicitly sets preconditioner)
89 -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
90 but use matrix-free approx for Jacobian-vector
91 products within Newton-Krylov method
92 */
93
94 PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
95
96 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 Customize nonlinear solver; set runtime options
98 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99
100 /*
101 Set an optional user-defined monitoring routine
102 */
103 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer));
104 PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0));
105
106 /*
107 Set names for some vectors to facilitate monitoring (optional)
108 */
109 PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
110 PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));
111
112 /*
113 Set SNES/KSP/KSP/PC runtime options, e.g.,
114 -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
115 */
116 PetscCall(SNESSetFromOptions(snes));
117
118 /*
119 Print parameters used for convergence testing (optional) ... just
120 to demonstrate this routine; this information is also printed with
121 the option -snes_view
122 */
123 PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
124 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));
125
126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 Initialize application:
128 Store right-hand side of PDE and exact solution
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130
131 xp = 0.0;
132 for (i = 0; i < n; i++) {
133 v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
134 PetscCall(VecSetValues(F, 1, &i, &v, INSERT_VALUES));
135 v = xp * xp * xp;
136 PetscCall(VecSetValues(U, 1, &i, &v, INSERT_VALUES));
137 xp += h;
138 }
139
140 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141 Evaluate initial guess; then solve nonlinear system
142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143 /*
144 Note: The user should initialize the vector, x, with the initial guess
145 for the nonlinear solver prior to calling SNESSolve(). In particular,
146 to employ an initial guess of zero, the user should explicitly set
147 this vector to zero by calling VecSet().
148 */
149 PetscCall(FormInitialGuess(x));
150 PetscCall(SNESSolve(snes, NULL, x));
151 PetscCall(SNESGetIterationNumber(snes, &its));
152 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its));
153
154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155 Check solution and clean up
156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157
158 /*
159 Check the error
160 */
161 PetscCall(VecAXPY(x, none, U));
162 PetscCall(VecNorm(x, NORM_2, &norm));
163 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
164
165 /*
166 Free work space. All PETSc objects should be destroyed when they
167 are no longer needed.
168 */
169 PetscCall(VecDestroy(&x));
170 PetscCall(VecDestroy(&r));
171 PetscCall(VecDestroy(&U));
172 PetscCall(VecDestroy(&F));
173 PetscCall(MatDestroy(&J));
174 PetscCall(SNESDestroy(&snes));
175 PetscCall(PetscViewerDestroy(&monP.viewer));
176 PetscCall(PetscFinalize());
177 return 0;
178 }
179 /* ------------------------------------------------------------------- */
180 /*
181 FormInitialGuess - Computes initial guess.
182
183 Input/Output Parameter:
184 . x - the solution vector
185 */
FormInitialGuess(Vec x)186 PetscErrorCode FormInitialGuess(Vec x)
187 {
188 PetscFunctionBeginUser;
189 PetscCall(VecSet(x, 0.5));
190 PetscFunctionReturn(PETSC_SUCCESS);
191 }
192 /* ------------------------------------------------------------------- */
193 /*
194 FormFunction - Evaluates nonlinear function, F(x).
195
196 Input Parameters:
197 . snes - the SNES context
198 . x - input vector
199 . ctx - optional user-defined context, as set by SNESSetFunction()
200
201 Output Parameter:
202 . f - function vector
203
204 Note:
205 The user-defined context can contain any application-specific data
206 needed for the function evaluation (such as various parameters, work
207 vectors, and grid information). In this program the context is just
208 a vector containing the right-hand side of the discretized PDE.
209 */
210
FormFunction(SNES snes,Vec x,Vec f,PetscCtx ctx)211 PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, PetscCtx ctx)
212 {
213 Vec g = (Vec)ctx;
214 const PetscScalar *xx, *gg;
215 PetscScalar *ff, d;
216 PetscInt i, n;
217
218 PetscFunctionBeginUser;
219 /*
220 Get pointers to vector data.
221 - For default PETSc vectors, VecGetArray() returns a pointer to
222 the data array. Otherwise, the routine is implementation dependent.
223 - You MUST call VecRestoreArray() when you no longer need access to
224 the array.
225 */
226 PetscCall(VecGetArrayRead(x, &xx));
227 PetscCall(VecGetArray(f, &ff));
228 PetscCall(VecGetArrayRead(g, &gg));
229
230 /*
231 Compute function
232 */
233 PetscCall(VecGetSize(x, &n));
234 d = (PetscReal)(n - 1);
235 d = d * d;
236 ff[0] = xx[0];
237 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];
238 ff[n - 1] = xx[n - 1] - 1.0;
239
240 /*
241 Restore vectors
242 */
243 PetscCall(VecRestoreArrayRead(x, &xx));
244 PetscCall(VecRestoreArray(f, &ff));
245 PetscCall(VecRestoreArrayRead(g, &gg));
246 PetscFunctionReturn(PETSC_SUCCESS);
247 }
248 /* ------------------------------------------------------------------- */
249 /*
250 FormJacobian - Evaluates Jacobian matrix.
251
252 Input Parameters:
253 . snes - the SNES context
254 . x - input vector
255 . dummy - optional user-defined context (not used here)
256
257 Output Parameters:
258 . jac - Jacobian matrix
259 . B - optionally different matrix used to construct the preconditioner
260
261 */
262
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void * dummy)263 PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
264 {
265 const PetscScalar *xx;
266 PetscScalar A[3], d;
267 PetscInt i, n, j[3];
268
269 PetscFunctionBeginUser;
270 /*
271 Get pointer to vector data
272 */
273 PetscCall(VecGetArrayRead(x, &xx));
274
275 /*
276 Compute Jacobian entries and insert into matrix.
277 - Note that in this case we set all elements for a particular
278 row at once.
279 */
280 PetscCall(VecGetSize(x, &n));
281 d = (PetscReal)(n - 1);
282 d = d * d;
283
284 /*
285 Interior grid points
286 */
287 for (i = 1; i < n - 1; i++) {
288 j[0] = i - 1;
289 j[1] = i;
290 j[2] = i + 1;
291 A[0] = A[2] = d;
292 A[1] = -2.0 * d + 2.0 * xx[i];
293 PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
294 }
295
296 /*
297 Boundary points
298 */
299 i = 0;
300 A[0] = 1.0;
301
302 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
303
304 i = n - 1;
305 A[0] = 1.0;
306
307 PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
308
309 /*
310 Restore vector
311 */
312 PetscCall(VecRestoreArrayRead(x, &xx));
313
314 /*
315 Assemble matrix
316 */
317 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
318 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
319 if (jac != B) {
320 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
321 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
322 }
323 PetscFunctionReturn(PETSC_SUCCESS);
324 }
325 /* ------------------------------------------------------------------- */
326 /*
327 Monitor - User-defined monitoring routine that views the
328 current iterate with an x-window plot.
329
330 Input Parameters:
331 snes - the SNES context
332 its - iteration number
333 norm - 2-norm function value (may be estimated)
334 ctx - optional user-defined context for private data for the
335 monitor routine, as set by SNESMonitorSet()
336
337 Note:
338 See the manpage for PetscViewerDrawOpen() for useful runtime options,
339 such as -nox to deactivate all x-window output.
340 */
Monitor(SNES snes,PetscInt its,PetscReal fnorm,PetscCtx ctx)341 PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, PetscCtx ctx)
342 {
343 MonitorCtx *monP = (MonitorCtx *)ctx;
344 Vec x;
345 SNESConvergedReason reason;
346
347 PetscFunctionBeginUser;
348 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm));
349 PetscCall(SNESGetConvergedReason(snes, &reason));
350 PetscCall(SNESGetSolution(snes, &x));
351 PetscCall(VecView(x, monP->viewer));
352 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " converged = %s\n", SNESConvergedReasons[reason]));
353 PetscFunctionReturn(PETSC_SUCCESS);
354 }
355
356 /*TEST
357
358 test:
359 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
360
361 test:
362 suffix: 2
363 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
364 requires: !single
365
366 test:
367 suffix: 3
368 args: -nox -malloc no -options_left no -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
369
370 test:
371 suffix: 4
372 args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
373 requires: !single
374
375 test:
376 suffix: 5
377 filter: grep -v atol | sed -e "s/CONVERGED_ITS/DIVERGED_MAX_IT/g" | sed -e "s/CONVERGED_FNORM_RELATIVE/DIVERGED_MAX_IT/g"
378 args: -nox -snes_type {{newtonls newtontr ncg ngmres qn anderson nrichardson ms ksponly ksptransposeonly vinewtonrsls vinewtonssls fas ms}} -snes_max_it 1
379 requires: !single
380
381 TEST*/
382