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: 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;
43: PetscFunctionBeginUser;
44: PetscCall(PetscInitialize(&argc, &argv, (char *)0, 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);
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Create nonlinear solver context
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
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));
69: /*
70: Set function evaluation routine and vector
71: */
72: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)F));
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create matrix data structure; set Jacobian evaluation routine
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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));
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 preconditioning matrix as set by the user,
90: but use matrix-free approx for Jacobian-vector
91: products within Newton-Krylov method
92: */
94: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Customize nonlinear solver; set runtime options
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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));
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"));
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));
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));
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Initialize application:
128: Store right-hand-side of PDE and exact solution
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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: }
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));
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Check solution and clean up
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
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));
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.
183: Input/Output Parameter:
184: . x - the solution vector
185: */
186: PetscErrorCode FormInitialGuess(Vec x)
187: {
188: PetscScalar pfive = .50;
189: PetscFunctionBeginUser;
190: PetscCall(VecSet(x, pfive));
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
193: /* ------------------------------------------------------------------- */
194: /*
195: FormFunction - Evaluates nonlinear function, F(x).
197: Input Parameters:
198: . snes - the SNES context
199: . x - input vector
200: . ctx - optional user-defined context, as set by SNESSetFunction()
202: Output Parameter:
203: . f - function vector
205: Note:
206: The user-defined context can contain any application-specific data
207: needed for the function evaluation (such as various parameters, work
208: vectors, and grid information). In this program the context is just
209: a vector containing the right-hand-side of the discretized PDE.
210: */
212: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
213: {
214: Vec g = (Vec)ctx;
215: const PetscScalar *xx, *gg;
216: PetscScalar *ff, d;
217: PetscInt i, n;
219: PetscFunctionBeginUser;
220: /*
221: Get pointers to vector data.
222: - For default PETSc vectors, VecGetArray() returns a pointer to
223: the data array. Otherwise, the routine is implementation dependent.
224: - You MUST call VecRestoreArray() when you no longer need access to
225: the array.
226: */
227: PetscCall(VecGetArrayRead(x, &xx));
228: PetscCall(VecGetArray(f, &ff));
229: PetscCall(VecGetArrayRead(g, &gg));
231: /*
232: Compute function
233: */
234: PetscCall(VecGetSize(x, &n));
235: d = (PetscReal)(n - 1);
236: d = d * d;
237: ff[0] = xx[0];
238: 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];
239: ff[n - 1] = xx[n - 1] - 1.0;
241: /*
242: Restore vectors
243: */
244: PetscCall(VecRestoreArrayRead(x, &xx));
245: PetscCall(VecRestoreArray(f, &ff));
246: PetscCall(VecRestoreArrayRead(g, &gg));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
249: /* ------------------------------------------------------------------- */
250: /*
251: FormJacobian - Evaluates Jacobian matrix.
253: Input Parameters:
254: . snes - the SNES context
255: . x - input vector
256: . dummy - optional user-defined context (not used here)
258: Output Parameters:
259: . jac - Jacobian matrix
260: . B - optionally different preconditioning matrix
262: */
264: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *dummy)
265: {
266: const PetscScalar *xx;
267: PetscScalar A[3], d;
268: PetscInt i, n, j[3];
270: PetscFunctionBeginUser;
271: /*
272: Get pointer to vector data
273: */
274: PetscCall(VecGetArrayRead(x, &xx));
276: /*
277: Compute Jacobian entries and insert into matrix.
278: - Note that in this case we set all elements for a particular
279: row at once.
280: */
281: PetscCall(VecGetSize(x, &n));
282: d = (PetscReal)(n - 1);
283: d = d * d;
285: /*
286: Interior grid points
287: */
288: for (i = 1; i < n - 1; i++) {
289: j[0] = i - 1;
290: j[1] = i;
291: j[2] = i + 1;
292: A[0] = A[2] = d;
293: A[1] = -2.0 * d + 2.0 * xx[i];
294: PetscCall(MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES));
295: }
297: /*
298: Boundary points
299: */
300: i = 0;
301: A[0] = 1.0;
303: PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
305: i = n - 1;
306: A[0] = 1.0;
308: PetscCall(MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES));
310: /*
311: Restore vector
312: */
313: PetscCall(VecRestoreArrayRead(x, &xx));
315: /*
316: Assemble matrix
317: */
318: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
319: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
320: if (jac != B) {
321: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
322: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
323: }
324: PetscFunctionReturn(PETSC_SUCCESS);
325: }
326: /* ------------------------------------------------------------------- */
327: /*
328: Monitor - User-defined monitoring routine that views the
329: current iterate with an x-window plot.
331: Input Parameters:
332: snes - the SNES context
333: its - iteration number
334: norm - 2-norm function value (may be estimated)
335: ctx - optional user-defined context for private data for the
336: monitor routine, as set by SNESMonitorSet()
338: Note:
339: See the manpage for PetscViewerDrawOpen() for useful runtime options,
340: such as -nox to deactivate all x-window output.
341: */
342: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
343: {
344: MonitorCtx *monP = (MonitorCtx *)ctx;
345: Vec x;
346: SNESConvergedReason reason;
348: PetscFunctionBeginUser;
349: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm));
350: PetscCall(SNESGetConvergedReason(snes, &reason));
351: PetscCall(SNESGetSolution(snes, &x));
352: PetscCall(VecView(x, monP->viewer));
353: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " converged = %s\n", SNESConvergedReasons[reason]));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*TEST
359: test:
360: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
362: test:
363: suffix: 2
364: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
365: requires: !single
367: test:
368: suffix: 3
369: 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
371: test:
372: suffix: 4
373: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
374: requires: !single
376: test:
377: suffix: 5
378: filter: grep -v atol | sed -e "s/CONVERGED_ITS/DIVERGED_MAX_IT/g" | sed -e "s/CONVERGED_FNORM_RELATIVE/DIVERGED_MAX_IT/g"
379: args: -nox -snes_type {{newtonls newtontr ncg ngmres qn anderson nrichardson ms ksponly ksptransposeonly vinewtonrsls vinewtonssls fas ms}} -snes_max_it 1
380: requires: !single
382: TEST*/