Actual source code: ex1.c
1: static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\
2: This example also illustrates the use of matrix coloring. Runtime options include:\n\
3: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
4: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
5: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
6: -my <yg>, where <yg> = number of grid points in the y-direction\n\n";
8: /*
10: Solid Fuel Ignition (SFI) problem. This problem is modeled by
11: the partial differential equation
13: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
15: with boundary conditions
17: u = 0 for x = 0, x = 1, y = 0, y = 1.
19: A finite difference approximation with the usual 5-point stencil
20: is used to discretize the boundary value problem to obtain a nonlinear
21: system of equations.
23: The parallel version of this code is snes/tutorials/ex5.c
25: */
27: /*
28: Include "petscsnes.h" so that we can use SNES solvers. Note that
29: this file automatically includes:
30: petscsys.h - base PETSc routines petscvec.h - vectors
31: petscmat.h - matrices
32: petscis.h - index sets petscksp.h - Krylov subspace methods
33: petscviewer.h - viewers petscpc.h - preconditioners
34: petscksp.h - linear solvers
35: */
37: #include <petscsnes.h>
39: /*
40: User-defined application context - contains data needed by the
41: application-provided call-back routines, FormJacobian() and
42: FormFunction().
43: */
44: typedef struct {
45: PetscReal param; /* test problem parameter */
46: PetscInt mx; /* Discretization in x-direction */
47: PetscInt my; /* Discretization in y-direction */
48: } AppCtx;
50: /*
51: User-defined routines
52: */
53: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
54: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
55: extern PetscErrorCode FormInitialGuess(AppCtx *, Vec);
56: extern PetscErrorCode ConvergenceTest(KSP, PetscInt, PetscReal, KSPConvergedReason *, void *);
57: extern PetscErrorCode ConvergenceDestroy(void *);
58: extern PetscErrorCode postcheck(SNES, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
60: int main(int argc, char **argv)
61: {
62: SNES snes; /* nonlinear solver context */
63: Vec x, r; /* solution, residual vectors */
64: Mat J; /* Jacobian matrix */
65: AppCtx user; /* user-defined application context */
66: PetscInt i, its, N, hist_its[50];
67: PetscMPIInt size;
68: PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0., history[50];
69: MatFDColoring fdcoloring;
70: PetscBool matrix_free = PETSC_FALSE, flg, fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE, pc = PETSC_FALSE, prunejacobian = PETSC_FALSE, null_appctx = PETSC_TRUE;
71: KSP ksp;
72: PetscInt *testarray;
74: PetscFunctionBeginUser;
75: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
76: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
77: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
79: /*
80: Initialize problem parameters
81: */
82: user.mx = 4;
83: user.my = 4;
84: user.param = 6.0;
85: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL));
86: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL));
87: PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL));
88: PetscCall(PetscOptionsGetBool(NULL, NULL, "-pc", &pc, NULL));
89: PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range");
90: N = user.mx * user.my;
91: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL));
92: PetscCall(PetscOptionsGetBool(NULL, NULL, "-prune_jacobian", &prunejacobian, NULL));
93: PetscCall(PetscOptionsGetBool(NULL, NULL, "-null_appctx", &null_appctx, NULL));
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create nonlinear solver context
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
101: if (pc) {
102: PetscCall(SNESSetType(snes, SNESNEWTONTR));
103: PetscCall(SNESNewtonTRSetPostCheck(snes, postcheck, NULL));
104: }
106: /* Test application context handling from Python */
107: if (!null_appctx) { PetscCall(SNESSetApplicationContext(snes, (void *)&user)); }
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create vector data structures; set function evaluation routine
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
114: PetscCall(VecSetSizes(x, PETSC_DECIDE, N));
115: PetscCall(VecSetFromOptions(x));
116: PetscCall(VecDuplicate(x, &r));
118: /*
119: Set function evaluation routine and vector. Whenever the nonlinear
120: solver needs to evaluate the nonlinear function, it will call this
121: routine.
122: - Note that the final routine argument is the user-defined
123: context that provides application-specific data for the
124: function evaluation routine.
125: */
126: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user));
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Create matrix data structure; set Jacobian evaluation routine
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: /*
133: Create matrix. Here we only approximately preallocate storage space
134: for the Jacobian. See the users manual for a discussion of better
135: techniques for preallocating matrix memory.
136: */
137: PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL));
138: if (!matrix_free) {
139: PetscBool matrix_free_operator = PETSC_FALSE;
140: PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &matrix_free_operator, NULL));
141: if (matrix_free_operator) matrix_free = PETSC_FALSE;
142: }
143: if (!matrix_free) PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, 5, NULL, &J));
145: /*
146: This option will cause the Jacobian to be computed via finite differences
147: efficiently using a coloring of the columns of the matrix.
148: */
149: PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_fd_coloring", &fd_coloring, NULL));
150: PetscCheck(!matrix_free || !fd_coloring, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Use only one of -snes_mf, -snes_fd_coloring options!\nYou can do -snes_mf_operator -snes_fd_coloring");
152: if (fd_coloring) {
153: ISColoring iscoloring;
154: MatColoring mc;
155: if (prunejacobian) {
156: /* Initialize x with random nonzero values so that the nonzeros in the Jacobian
157: can better reflect the sparsity structure of the Jacobian. */
158: PetscRandom rctx;
159: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
160: PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0));
161: PetscCall(VecSetRandom(x, rctx));
162: PetscCall(PetscRandomDestroy(&rctx));
163: }
165: /*
166: This initializes the nonzero structure of the Jacobian. This is artificial
167: because clearly if we had a routine to compute the Jacobian we won't need
168: to use finite differences.
169: */
170: PetscCall(FormJacobian(snes, x, J, J, &user));
172: /*
173: Color the matrix, i.e. determine groups of columns that share no common
174: rows. These columns in the Jacobian can all be computed simultaneously.
175: */
176: PetscCall(MatColoringCreate(J, &mc));
177: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
178: PetscCall(MatColoringSetFromOptions(mc));
179: PetscCall(MatColoringApply(mc, &iscoloring));
180: PetscCall(MatColoringDestroy(&mc));
181: /*
182: Create the data structure that SNESComputeJacobianDefaultColor() uses
183: to compute the actual Jacobians via finite differences.
184: */
185: PetscCall(MatFDColoringCreate(J, iscoloring, &fdcoloring));
186: PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))FormFunction, &user));
187: PetscCall(MatFDColoringSetFromOptions(fdcoloring));
188: PetscCall(MatFDColoringSetUp(J, iscoloring, fdcoloring));
189: /*
190: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
191: to compute Jacobians.
192: */
193: PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring));
194: PetscCall(ISColoringDestroy(&iscoloring));
195: if (prunejacobian) PetscCall(SNESPruneJacobianColor(snes, J, J));
196: }
197: /*
198: Set Jacobian matrix data structure and default Jacobian evaluation
199: routine. Whenever the nonlinear solver needs to compute the
200: Jacobian matrix, it will call this routine.
201: - Note that the final routine argument is the user-defined
202: context that provides application-specific data for the
203: Jacobian evaluation routine.
204: - The user can override with:
205: -snes_fd : default finite differencing approximation of Jacobian
206: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
207: (unless user explicitly sets preconditioner)
208: -snes_mf_operator : form preconditioning matrix as set by the user,
209: but use matrix-free approx for Jacobian-vector
210: products within Newton-Krylov method
211: */
212: else if (!matrix_free) {
213: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, (void *)&user));
214: }
216: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: Customize nonlinear solver; set runtime options
218: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
220: /*
221: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
222: */
223: PetscCall(SNESSetFromOptions(snes));
225: /*
226: Set array that saves the function norms. This array is intended
227: when the user wants to save the convergence history for later use
228: rather than just to view the function norms via -snes_monitor.
229: */
230: PetscCall(SNESSetConvergenceHistory(snes, history, hist_its, 50, PETSC_TRUE));
232: /*
233: Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
234: user provided test before the specialized test. The convergence context is just an array to
235: test that it gets properly freed at the end
236: */
237: if (use_convergence_test) {
238: PetscCall(SNESGetKSP(snes, &ksp));
239: PetscCall(PetscMalloc1(5, &testarray));
240: PetscCall(KSPSetConvergenceTest(ksp, ConvergenceTest, testarray, ConvergenceDestroy));
241: }
243: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244: Evaluate initial guess; then solve nonlinear system
245: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
246: /*
247: Note: The user should initialize the vector, x, with the initial guess
248: for the nonlinear solver prior to calling SNESSolve(). In particular,
249: to employ an initial guess of zero, the user should explicitly set
250: this vector to zero by calling VecSet().
251: */
252: PetscCall(FormInitialGuess(&user, x));
253: PetscCall(SNESSolve(snes, NULL, x));
254: PetscCall(SNESGetIterationNumber(snes, &its));
255: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
257: /*
258: Print the convergence history. This is intended just to demonstrate
259: use of the data attained via SNESSetConvergenceHistory().
260: */
261: PetscCall(PetscOptionsHasName(NULL, NULL, "-print_history", &flg));
262: if (flg) {
263: for (i = 0; i < its + 1; i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iteration %" PetscInt_FMT ": Linear iterations %" PetscInt_FMT " Function norm = %g\n", i, hist_its[i], (double)history[i]));
264: }
266: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267: Free work space. All PETSc objects should be destroyed when they
268: are no longer needed.
269: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: if (!matrix_free) PetscCall(MatDestroy(&J));
272: if (fd_coloring) PetscCall(MatFDColoringDestroy(&fdcoloring));
273: PetscCall(VecDestroy(&x));
274: PetscCall(VecDestroy(&r));
275: PetscCall(SNESDestroy(&snes));
276: PetscCall(PetscFinalize());
277: return 0;
278: }
280: /*
281: FormInitialGuess - Forms initial approximation.
283: Input Parameters:
284: user - user-defined application context
285: X - vector
287: Output Parameter:
288: X - vector
289: */
290: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
291: {
292: PetscInt i, j, row, mx, my;
293: PetscReal lambda, temp1, temp, hx, hy;
294: PetscScalar *x;
296: PetscFunctionBeginUser;
297: mx = user->mx;
298: my = user->my;
299: lambda = user->param;
301: hx = 1.0 / (PetscReal)(mx - 1);
302: hy = 1.0 / (PetscReal)(my - 1);
304: /*
305: Get a pointer to vector data.
306: - For default PETSc vectors, VecGetArray() returns a pointer to
307: the data array. Otherwise, the routine is implementation dependent.
308: - You MUST call VecRestoreArray() when you no longer need access to
309: the array.
310: */
311: PetscCall(VecGetArray(X, &x));
312: temp1 = lambda / (lambda + 1.0);
313: for (j = 0; j < my; j++) {
314: temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
315: for (i = 0; i < mx; i++) {
316: row = i + j * mx;
317: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
318: x[row] = 0.0;
319: continue;
320: }
321: x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
322: }
323: }
325: /*
326: Restore vector
327: */
328: PetscCall(VecRestoreArray(X, &x));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /*
333: FormFunction - Evaluates nonlinear function, F(x).
335: Input Parameters:
336: . snes - the SNES context
337: . X - input vector
338: . ptr - optional user-defined context, as set by SNESSetFunction()
340: Output Parameter:
341: . F - function vector
342: */
343: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
344: {
345: AppCtx *user = (AppCtx *)ptr;
346: PetscInt i, j, row, mx, my;
347: PetscReal two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx;
348: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f;
349: const PetscScalar *x;
351: PetscFunctionBeginUser;
352: mx = user->mx;
353: my = user->my;
354: lambda = user->param;
355: hx = one / (PetscReal)(mx - 1);
356: hy = one / (PetscReal)(my - 1);
357: sc = hx * hy;
358: hxdhy = hx / hy;
359: hydhx = hy / hx;
361: /*
362: Get pointers to vector data
363: */
364: PetscCall(VecGetArrayRead(X, &x));
365: PetscCall(VecGetArray(F, &f));
367: /*
368: Compute function
369: */
370: for (j = 0; j < my; j++) {
371: for (i = 0; i < mx; i++) {
372: row = i + j * mx;
373: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
374: f[row] = x[row];
375: continue;
376: }
377: u = x[row];
378: ub = x[row - mx];
379: ul = x[row - 1];
380: ut = x[row + mx];
381: ur = x[row + 1];
382: uxx = (-ur + two * u - ul) * hydhx;
383: uyy = (-ut + two * u - ub) * hxdhy;
384: f[row] = uxx + uyy - sc * lambda * PetscExpScalar(u);
385: }
386: }
388: /*
389: Restore vectors
390: */
391: PetscCall(VecRestoreArrayRead(X, &x));
392: PetscCall(VecRestoreArray(F, &f));
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*
397: FormJacobian - Evaluates Jacobian matrix.
399: Input Parameters:
400: . snes - the SNES context
401: . x - input vector
402: . ptr - optional user-defined context, as set by SNESSetJacobian()
404: Output Parameters:
405: . A - Jacobian matrix
406: . B - optionally different preconditioning matrix
407: . flag - flag indicating matrix structure
408: */
409: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
410: {
411: AppCtx *user = (AppCtx *)ptr; /* user-defined application context */
412: PetscInt i, j, row, mx, my, col[5];
413: PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc;
414: const PetscScalar *x;
415: PetscReal hx, hy, hxdhy, hydhx;
417: PetscFunctionBeginUser;
418: mx = user->mx;
419: my = user->my;
420: lambda = user->param;
421: hx = 1.0 / (PetscReal)(mx - 1);
422: hy = 1.0 / (PetscReal)(my - 1);
423: sc = hx * hy;
424: hxdhy = hx / hy;
425: hydhx = hy / hx;
427: /*
428: Get pointer to vector data
429: */
430: PetscCall(VecGetArrayRead(X, &x));
432: /*
433: Compute entries of the Jacobian
434: */
435: for (j = 0; j < my; j++) {
436: for (i = 0; i < mx; i++) {
437: row = i + j * mx;
438: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
439: PetscCall(MatSetValues(jac, 1, &row, 1, &row, &one, INSERT_VALUES));
440: continue;
441: }
442: v[0] = -hxdhy;
443: col[0] = row - mx;
444: v[1] = -hydhx;
445: col[1] = row - 1;
446: v[2] = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]);
447: col[2] = row;
448: v[3] = -hydhx;
449: col[3] = row + 1;
450: v[4] = -hxdhy;
451: col[4] = row + mx;
452: PetscCall(MatSetValues(jac, 1, &row, 5, col, v, INSERT_VALUES));
453: }
454: }
456: /*
457: Restore vector
458: */
459: PetscCall(VecRestoreArrayRead(X, &x));
461: /*
462: Assemble matrix
463: */
464: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
465: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
467: if (jac != J) {
468: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
469: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
470: }
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: PetscErrorCode ConvergenceTest(KSP ksp, PetscInt it, PetscReal nrm, KSPConvergedReason *reason, void *ctx)
476: {
477: PetscFunctionBeginUser;
478: *reason = KSP_CONVERGED_ITERATING;
479: if (it > 1) {
480: *reason = KSP_CONVERGED_ITS;
481: PetscCall(PetscInfo(NULL, "User provided convergence test returning after 2 iterations\n"));
482: }
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: PetscErrorCode ConvergenceDestroy(void *ctx)
487: {
488: PetscFunctionBeginUser;
489: PetscCall(PetscInfo(NULL, "User provided convergence destroy called\n"));
490: PetscCall(PetscFree(ctx));
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: PetscErrorCode postcheck(SNES snes, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, void *ctx)
495: {
496: PetscReal norm;
497: Vec tmp;
499: PetscFunctionBeginUser;
500: PetscCall(VecDuplicate(x, &tmp));
501: PetscCall(VecWAXPY(tmp, -1.0, x, w));
502: PetscCall(VecNorm(tmp, NORM_2, &norm));
503: PetscCall(VecDestroy(&tmp));
504: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of search step %g\n", (double)norm));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: /*TEST
510: build:
511: requires: !single
513: test:
514: args: -ksp_gmres_cgs_refinement_type refine_always
516: test:
517: suffix: 2
518: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
520: test:
521: suffix: 2a
522: filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
523: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
524: requires: defined(PETSC_USE_INFO)
526: test:
527: suffix: 2b
528: filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test"
529: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
530: requires: defined(PETSC_USE_INFO)
532: test:
533: suffix: 3
534: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
536: test:
537: suffix: 4
538: args: -pc -par 6.807 -snes_monitor -snes_converged_reason
540: test:
541: suffix: 5
542: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always -prune_jacobian
543: output_file: output/ex1_3.out
545: test:
546: suffix: 6
547: args: -snes_monitor draw:image:testfile -viewer_view
549: test:
550: suffix: python
551: requires: petsc4py
552: args: -python -snes_type python -snes_python_type ex1.py:MySNES -snes_view -null_appctx {{0 1}separate output}
553: localrunfiles: ex1.py
555: TEST*/