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;
72: KSP ksp;
73: PetscInt *testarray;
76: PetscInitialize(&argc, &argv, (char *)0, help);
77: MPI_Comm_size(PETSC_COMM_WORLD, &size);
80: /*
81: Initialize problem parameters
82: */
83: user.mx = 4;
84: user.my = 4;
85: user.param = 6.0;
86: PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL);
87: PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL);
88: PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL);
89: PetscOptionsGetBool(NULL, NULL, "-pc", &pc, NULL);
91: N = user.mx * user.my;
92: PetscOptionsGetBool(NULL, NULL, "-use_convergence_test", &use_convergence_test, NULL);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create nonlinear solver context
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: SNESCreate(PETSC_COMM_WORLD, &snes);
100: if (pc) {
101: SNESSetType(snes, SNESNEWTONTR);
102: SNESNewtonTRSetPostCheck(snes, postcheck, NULL);
103: }
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create vector data structures; set function evaluation routine
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: VecCreate(PETSC_COMM_WORLD, &x);
110: VecSetSizes(x, PETSC_DECIDE, N);
111: VecSetFromOptions(x);
112: VecDuplicate(x, &r);
114: /*
115: Set function evaluation routine and vector. Whenever the nonlinear
116: solver needs to evaluate the nonlinear function, it will call this
117: routine.
118: - Note that the final routine argument is the user-defined
119: context that provides application-specific data for the
120: function evaluation routine.
121: */
122: SNESSetFunction(snes, r, FormFunction, (void *)&user);
124: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Create matrix data structure; set Jacobian evaluation routine
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: /*
129: Create matrix. Here we only approximately preallocate storage space
130: for the Jacobian. See the users manual for a discussion of better
131: techniques for preallocating matrix memory.
132: */
133: PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL);
134: if (!matrix_free) {
135: PetscBool matrix_free_operator = PETSC_FALSE;
136: PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &matrix_free_operator, NULL);
137: if (matrix_free_operator) matrix_free = PETSC_FALSE;
138: }
139: if (!matrix_free) MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, 5, NULL, &J);
141: /*
142: This option will cause the Jacobian to be computed via finite differences
143: efficiently using a coloring of the columns of the matrix.
144: */
145: PetscOptionsGetBool(NULL, NULL, "-snes_fd_coloring", &fd_coloring, NULL);
148: if (fd_coloring) {
149: ISColoring iscoloring;
150: MatColoring mc;
152: /*
153: This initializes the nonzero structure of the Jacobian. This is artificial
154: because clearly if we had a routine to compute the Jacobian we won't need
155: to use finite differences.
156: */
157: FormJacobian(snes, x, J, J, &user);
159: /*
160: Color the matrix, i.e. determine groups of columns that share no common
161: rows. These columns in the Jacobian can all be computed simultaneously.
162: */
163: MatColoringCreate(J, &mc);
164: MatColoringSetType(mc, MATCOLORINGSL);
165: MatColoringSetFromOptions(mc);
166: MatColoringApply(mc, &iscoloring);
167: MatColoringDestroy(&mc);
168: /*
169: Create the data structure that SNESComputeJacobianDefaultColor() uses
170: to compute the actual Jacobians via finite differences.
171: */
172: MatFDColoringCreate(J, iscoloring, &fdcoloring);
173: MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))FormFunction, &user);
174: MatFDColoringSetFromOptions(fdcoloring);
175: MatFDColoringSetUp(J, iscoloring, fdcoloring);
176: /*
177: Tell SNES to use the routine SNESComputeJacobianDefaultColor()
178: to compute Jacobians.
179: */
180: SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring);
181: ISColoringDestroy(&iscoloring);
182: }
183: /*
184: Set Jacobian matrix data structure and default Jacobian evaluation
185: routine. Whenever the nonlinear solver needs to compute the
186: Jacobian matrix, it will call this routine.
187: - Note that the final routine argument is the user-defined
188: context that provides application-specific data for the
189: Jacobian evaluation routine.
190: - The user can override with:
191: -snes_fd : default finite differencing approximation of Jacobian
192: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
193: (unless user explicitly sets preconditioner)
194: -snes_mf_operator : form preconditioning matrix as set by the user,
195: but use matrix-free approx for Jacobian-vector
196: products within Newton-Krylov method
197: */
198: else if (!matrix_free) {
199: SNESSetJacobian(snes, J, J, FormJacobian, (void *)&user);
200: }
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Customize nonlinear solver; set runtime options
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: /*
207: Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
208: */
209: SNESSetFromOptions(snes);
211: /*
212: Set array that saves the function norms. This array is intended
213: when the user wants to save the convergence history for later use
214: rather than just to view the function norms via -snes_monitor.
215: */
216: SNESSetConvergenceHistory(snes, history, hist_its, 50, PETSC_TRUE);
218: /*
219: Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
220: user provided test before the specialized test. The convergence context is just an array to
221: test that it gets properly freed at the end
222: */
223: if (use_convergence_test) {
224: SNESGetKSP(snes, &ksp);
225: PetscMalloc1(5, &testarray);
226: KSPSetConvergenceTest(ksp, ConvergenceTest, testarray, ConvergenceDestroy);
227: }
229: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230: Evaluate initial guess; then solve nonlinear system
231: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232: /*
233: Note: The user should initialize the vector, x, with the initial guess
234: for the nonlinear solver prior to calling SNESSolve(). In particular,
235: to employ an initial guess of zero, the user should explicitly set
236: this vector to zero by calling VecSet().
237: */
238: FormInitialGuess(&user, x);
239: SNESSolve(snes, NULL, x);
240: SNESGetIterationNumber(snes, &its);
241: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
243: /*
244: Print the convergence history. This is intended just to demonstrate
245: use of the data attained via SNESSetConvergenceHistory().
246: */
247: PetscOptionsHasName(NULL, NULL, "-print_history", &flg);
248: if (flg) {
249: for (i = 0; i < its + 1; i++) PetscPrintf(PETSC_COMM_WORLD, "iteration %" PetscInt_FMT ": Linear iterations %" PetscInt_FMT " Function norm = %g\n", i, hist_its[i], (double)history[i]);
250: }
252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: Free work space. All PETSc objects should be destroyed when they
254: are no longer needed.
255: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
257: if (!matrix_free) MatDestroy(&J);
258: if (fd_coloring) MatFDColoringDestroy(&fdcoloring);
259: VecDestroy(&x);
260: VecDestroy(&r);
261: SNESDestroy(&snes);
262: PetscFinalize();
263: return 0;
264: }
266: /*
267: FormInitialGuess - Forms initial approximation.
269: Input Parameters:
270: user - user-defined application context
271: X - vector
273: Output Parameter:
274: X - vector
275: */
276: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
277: {
278: PetscInt i, j, row, mx, my;
279: PetscReal lambda, temp1, temp, hx, hy;
280: PetscScalar *x;
282: mx = user->mx;
283: my = user->my;
284: lambda = user->param;
286: hx = 1.0 / (PetscReal)(mx - 1);
287: hy = 1.0 / (PetscReal)(my - 1);
289: /*
290: Get a pointer to vector data.
291: - For default PETSc vectors, VecGetArray() returns a pointer to
292: the data array. Otherwise, the routine is implementation dependent.
293: - You MUST call VecRestoreArray() when you no longer need access to
294: the array.
295: */
296: VecGetArray(X, &x);
297: temp1 = lambda / (lambda + 1.0);
298: for (j = 0; j < my; j++) {
299: temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
300: for (i = 0; i < mx; i++) {
301: row = i + j * mx;
302: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
303: x[row] = 0.0;
304: continue;
305: }
306: x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
307: }
308: }
310: /*
311: Restore vector
312: */
313: VecRestoreArray(X, &x);
314: return 0;
315: }
317: /*
318: FormFunction - Evaluates nonlinear function, F(x).
320: Input Parameters:
321: . snes - the SNES context
322: . X - input vector
323: . ptr - optional user-defined context, as set by SNESSetFunction()
325: Output Parameter:
326: . F - function vector
327: */
328: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
329: {
330: AppCtx *user = (AppCtx *)ptr;
331: PetscInt i, j, row, mx, my;
332: PetscReal two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx;
333: PetscScalar ut, ub, ul, ur, u, uxx, uyy, sc, *f;
334: const PetscScalar *x;
336: mx = user->mx;
337: my = user->my;
338: lambda = user->param;
339: hx = one / (PetscReal)(mx - 1);
340: hy = one / (PetscReal)(my - 1);
341: sc = hx * hy;
342: hxdhy = hx / hy;
343: hydhx = hy / hx;
345: /*
346: Get pointers to vector data
347: */
348: VecGetArrayRead(X, &x);
349: VecGetArray(F, &f);
351: /*
352: Compute function
353: */
354: for (j = 0; j < my; j++) {
355: for (i = 0; i < mx; i++) {
356: row = i + j * mx;
357: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
358: f[row] = x[row];
359: continue;
360: }
361: u = x[row];
362: ub = x[row - mx];
363: ul = x[row - 1];
364: ut = x[row + mx];
365: ur = x[row + 1];
366: uxx = (-ur + two * u - ul) * hydhx;
367: uyy = (-ut + two * u - ub) * hxdhy;
368: f[row] = uxx + uyy - sc * lambda * PetscExpScalar(u);
369: }
370: }
372: /*
373: Restore vectors
374: */
375: VecRestoreArrayRead(X, &x);
376: VecRestoreArray(F, &f);
377: return 0;
378: }
380: /*
381: FormJacobian - Evaluates Jacobian matrix.
383: Input Parameters:
384: . snes - the SNES context
385: . x - input vector
386: . ptr - optional user-defined context, as set by SNESSetJacobian()
388: Output Parameters:
389: . A - Jacobian matrix
390: . B - optionally different preconditioning matrix
391: . flag - flag indicating matrix structure
392: */
393: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
394: {
395: AppCtx *user = (AppCtx *)ptr; /* user-defined application context */
396: PetscInt i, j, row, mx, my, col[5];
397: PetscScalar two = 2.0, one = 1.0, lambda, v[5], sc;
398: const PetscScalar *x;
399: PetscReal hx, hy, hxdhy, hydhx;
401: mx = user->mx;
402: my = user->my;
403: lambda = user->param;
404: hx = 1.0 / (PetscReal)(mx - 1);
405: hy = 1.0 / (PetscReal)(my - 1);
406: sc = hx * hy;
407: hxdhy = hx / hy;
408: hydhx = hy / hx;
410: /*
411: Get pointer to vector data
412: */
413: VecGetArrayRead(X, &x);
415: /*
416: Compute entries of the Jacobian
417: */
418: for (j = 0; j < my; j++) {
419: for (i = 0; i < mx; i++) {
420: row = i + j * mx;
421: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
422: MatSetValues(jac, 1, &row, 1, &row, &one, INSERT_VALUES);
423: continue;
424: }
425: v[0] = -hxdhy;
426: col[0] = row - mx;
427: v[1] = -hydhx;
428: col[1] = row - 1;
429: v[2] = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]);
430: col[2] = row;
431: v[3] = -hydhx;
432: col[3] = row + 1;
433: v[4] = -hxdhy;
434: col[4] = row + mx;
435: MatSetValues(jac, 1, &row, 5, col, v, INSERT_VALUES);
436: }
437: }
439: /*
440: Restore vector
441: */
442: VecRestoreArrayRead(X, &x);
444: /*
445: Assemble matrix
446: */
447: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
448: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
450: if (jac != J) {
451: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
452: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
453: }
455: return 0;
456: }
458: PetscErrorCode ConvergenceTest(KSP ksp, PetscInt it, PetscReal nrm, KSPConvergedReason *reason, void *ctx)
459: {
460: *reason = KSP_CONVERGED_ITERATING;
461: if (it > 1) {
462: *reason = KSP_CONVERGED_ITS;
463: PetscInfo(NULL, "User provided convergence test returning after 2 iterations\n");
464: }
465: return 0;
466: }
468: PetscErrorCode ConvergenceDestroy(void *ctx)
469: {
470: PetscInfo(NULL, "User provided convergence destroy called\n");
471: PetscFree(ctx);
472: return 0;
473: }
475: PetscErrorCode postcheck(SNES snes, Vec x, Vec y, Vec w, PetscBool *changed_y, PetscBool *changed_w, void *ctx)
476: {
477: PetscReal norm;
478: Vec tmp;
480: VecDuplicate(x, &tmp);
481: VecWAXPY(tmp, -1.0, x, w);
482: VecNorm(tmp, NORM_2, &norm);
483: VecDestroy(&tmp);
484: PetscPrintf(PETSC_COMM_WORLD, "Norm of search step %g\n", (double)norm);
485: return 0;
486: }
488: /*TEST
490: build:
491: requires: !single
493: test:
494: args: -ksp_gmres_cgs_refinement_type refine_always
496: test:
497: suffix: 2
498: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always
500: test:
501: suffix: 2a
502: filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
503: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
504: requires: defined(PETSC_USE_INFO)
506: test:
507: suffix: 2b
508: filter: grep -i "User provided convergence test" > /dev/null && echo "Found User provided convergence test"
509: args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
510: requires: defined(PETSC_USE_INFO)
512: test:
513: suffix: 3
514: args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always
516: test:
517: suffix: 4
518: args: -pc -par 6.807 -snes_monitor -snes_converged_reason
520: TEST*/