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