Actual source code: ex2.c
2: static char help[] = "Newton method to solve u'' + u^{2} = f, sequentially.\n\
3: This example employs a user-defined monitoring routine.\n\n";
5: /*
6: Include "petscdraw.h" so that we can use PETSc drawing routines.
7: Include "petscsnes.h" so that we can use SNES solvers. Note that this
8: file automatically includes:
9: petscsys.h - base PETSc routines petscvec.h - vectors
10: petscmat.h - matrices
11: petscis.h - index sets petscksp.h - Krylov subspace methods
12: petscviewer.h - viewers petscpc.h - preconditioners
13: petscksp.h - linear solvers
14: */
16: #include <petscsnes.h>
18: /*
19: User-defined routines
20: */
21: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
22: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
23: extern PetscErrorCode FormInitialGuess(Vec);
24: extern PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
26: /*
27: User-defined context for monitoring
28: */
29: typedef struct {
30: PetscViewer viewer;
31: } MonitorCtx;
33: int main(int argc, char **argv)
34: {
35: SNES snes; /* SNES context */
36: Vec x, r, F, U; /* vectors */
37: Mat J; /* Jacobian matrix */
38: MonitorCtx monP; /* monitoring context */
39: PetscInt its, n = 5, i, maxit, maxf;
40: PetscMPIInt size;
41: PetscScalar h, xp, v, none = -1.0;
42: PetscReal abstol, rtol, stol, norm;
45: PetscInitialize(&argc, &argv, (char *)0, help);
46: MPI_Comm_size(PETSC_COMM_WORLD, &size);
48: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
49: h = 1.0 / (n - 1);
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create nonlinear solver context
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: SNESCreate(PETSC_COMM_WORLD, &snes);
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create vector data structures; set function evaluation routine
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: /*
61: Note that we form 1 vector from scratch and then duplicate as needed.
62: */
63: VecCreate(PETSC_COMM_WORLD, &x);
64: VecSetSizes(x, PETSC_DECIDE, n);
65: VecSetFromOptions(x);
66: VecDuplicate(x, &r);
67: VecDuplicate(x, &F);
68: VecDuplicate(x, &U);
70: /*
71: Set function evaluation routine and vector
72: */
73: SNESSetFunction(snes, r, FormFunction, (void *)F);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create matrix data structure; set Jacobian evaluation routine
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: MatCreate(PETSC_COMM_WORLD, &J);
80: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, n, n);
81: MatSetFromOptions(J);
82: MatSeqAIJSetPreallocation(J, 3, NULL);
84: /*
85: Set Jacobian matrix data structure and default Jacobian evaluation
86: routine. User can override with:
87: -snes_fd : default finite differencing approximation of Jacobian
88: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
89: (unless user explicitly sets preconditioner)
90: -snes_mf_operator : form preconditioning matrix as set by the user,
91: but use matrix-free approx for Jacobian-vector
92: products within Newton-Krylov method
93: */
95: SNESSetJacobian(snes, J, J, FormJacobian, NULL);
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Customize nonlinear solver; set runtime options
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: /*
102: Set an optional user-defined monitoring routine
103: */
104: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer);
105: SNESMonitorSet(snes, Monitor, &monP, 0);
107: /*
108: Set names for some vectors to facilitate monitoring (optional)
109: */
110: PetscObjectSetName((PetscObject)x, "Approximate Solution");
111: PetscObjectSetName((PetscObject)U, "Exact Solution");
113: /*
114: Set SNES/KSP/KSP/PC runtime options, e.g.,
115: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
116: */
117: SNESSetFromOptions(snes);
119: /*
120: Print parameters used for convergence testing (optional) ... just
121: to demonstrate this routine; this information is also printed with
122: the option -snes_view
123: */
124: SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf);
125: 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);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Initialize application:
129: Store right-hand-side of PDE and exact solution
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: xp = 0.0;
133: for (i = 0; i < n; i++) {
134: v = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
135: VecSetValues(F, 1, &i, &v, INSERT_VALUES);
136: v = xp * xp * xp;
137: VecSetValues(U, 1, &i, &v, INSERT_VALUES);
138: xp += h;
139: }
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Evaluate initial guess; then solve nonlinear system
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: /*
145: Note: The user should initialize the vector, x, with the initial guess
146: for the nonlinear solver prior to calling SNESSolve(). In particular,
147: to employ an initial guess of zero, the user should explicitly set
148: this vector to zero by calling VecSet().
149: */
150: FormInitialGuess(x);
151: SNESSolve(snes, NULL, x);
152: SNESGetIterationNumber(snes, &its);
153: PetscPrintf(PETSC_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT "\n\n", its);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Check solution and clean up
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: /*
160: Check the error
161: */
162: VecAXPY(x, none, U);
163: VecNorm(x, NORM_2, &norm);
164: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
166: /*
167: Free work space. All PETSc objects should be destroyed when they
168: are no longer needed.
169: */
170: VecDestroy(&x);
171: VecDestroy(&r);
172: VecDestroy(&U);
173: VecDestroy(&F);
174: MatDestroy(&J);
175: SNESDestroy(&snes);
176: PetscViewerDestroy(&monP.viewer);
177: PetscFinalize();
178: return 0;
179: }
180: /* ------------------------------------------------------------------- */
181: /*
182: FormInitialGuess - Computes initial guess.
184: Input/Output Parameter:
185: . x - the solution vector
186: */
187: PetscErrorCode FormInitialGuess(Vec x)
188: {
189: PetscScalar pfive = .50;
190: VecSet(x, pfive);
191: return 0;
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: /*
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: VecGetArrayRead(x, &xx);
227: VecGetArray(f, &ff);
228: VecGetArrayRead(g, &gg);
230: /*
231: Compute function
232: */
233: 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;
240: /*
241: Restore vectors
242: */
243: VecRestoreArrayRead(x, &xx);
244: VecRestoreArray(f, &ff);
245: VecRestoreArrayRead(g, &gg);
246: return 0;
247: }
248: /* ------------------------------------------------------------------- */
249: /*
250: FormJacobian - Evaluates Jacobian matrix.
252: Input Parameters:
253: . snes - the SNES context
254: . x - input vector
255: . dummy - optional user-defined context (not used here)
257: Output Parameters:
258: . jac - Jacobian matrix
259: . B - optionally different preconditioning matrix
261: */
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];
269: /*
270: Get pointer to vector data
271: */
272: VecGetArrayRead(x, &xx);
274: /*
275: Compute Jacobian entries and insert into matrix.
276: - Note that in this case we set all elements for a particular
277: row at once.
278: */
279: VecGetSize(x, &n);
280: d = (PetscReal)(n - 1);
281: d = d * d;
283: /*
284: Interior grid points
285: */
286: for (i = 1; i < n - 1; i++) {
287: j[0] = i - 1;
288: j[1] = i;
289: j[2] = i + 1;
290: A[0] = A[2] = d;
291: A[1] = -2.0 * d + 2.0 * xx[i];
292: MatSetValues(B, 1, &i, 3, j, A, INSERT_VALUES);
293: }
295: /*
296: Boundary points
297: */
298: i = 0;
299: A[0] = 1.0;
301: MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES);
303: i = n - 1;
304: A[0] = 1.0;
306: MatSetValues(B, 1, &i, 1, &i, A, INSERT_VALUES);
308: /*
309: Restore vector
310: */
311: VecRestoreArrayRead(x, &xx);
313: /*
314: Assemble matrix
315: */
316: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
317: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
318: if (jac != B) {
319: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
320: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
321: }
322: return 0;
323: }
324: /* ------------------------------------------------------------------- */
325: /*
326: Monitor - User-defined monitoring routine that views the
327: current iterate with an x-window plot.
329: Input Parameters:
330: snes - the SNES context
331: its - iteration number
332: norm - 2-norm function value (may be estimated)
333: ctx - optional user-defined context for private data for the
334: monitor routine, as set by SNESMonitorSet()
336: Note:
337: See the manpage for PetscViewerDrawOpen() for useful runtime options,
338: such as -nox to deactivate all x-window output.
339: */
340: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
341: {
342: MonitorCtx *monP = (MonitorCtx *)ctx;
343: Vec x;
345: PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ", SNES Function norm %g\n", its, (double)fnorm);
346: SNESGetSolution(snes, &x);
347: VecView(x, monP->viewer);
348: return 0;
349: }
351: /*TEST
353: test:
354: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_view -pc_type jacobi -ksp_gmres_cgs_refinement_type refine_always
356: test:
357: suffix: 2
358: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontr -snes_view
359: requires: !single
361: test:
362: suffix: 3
363: 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
365: test:
366: suffix: 4
367: args: -nox -snes_monitor_cancel -snes_monitor_short -snes_type newtontrdc -snes_view
368: requires: !single
370: TEST*/