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