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