Actual source code: ex14.c
2: static char help[] = "Bratu nonlinear PDE in 3d.\n\
3: We solve the Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\
4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
5: The command line options include:\n\
6: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
7: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\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, z = 0, z = 1
20: A finite difference approximation with the usual 7-point stencil
21: is used to discretize the boundary value problem to obtain a nonlinear
22: system of equations.
24: ------------------------------------------------------------------------- */
26: /*
27: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
28: Include "petscsnes.h" so that we can use SNES solvers. Note that this
29: 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: */
36: #include <petscdm.h>
37: #include <petscdmda.h>
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: DM da; /* distributed array data structure */
48: } AppCtx;
50: /*
51: User-defined routines
52: */
53: extern PetscErrorCode FormFunctionLocal(SNES, Vec, Vec, void *);
54: extern PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
55: extern PetscErrorCode FormInitialGuess(AppCtx *, Vec);
56: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
58: int main(int argc, char **argv)
59: {
60: SNES snes; /* nonlinear solver */
61: Vec x, r; /* solution, residual vectors */
62: Mat J = NULL; /* Jacobian matrix */
63: AppCtx user; /* user-defined work context */
64: PetscInt its; /* iterations for convergence */
65: MatFDColoring matfdcoloring = NULL;
66: PetscBool matrix_free = PETSC_FALSE, coloring = PETSC_FALSE, coloring_ds = PETSC_FALSE, local_coloring = PETSC_FALSE;
67: PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0., fnorm;
69: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70: Initialize program
71: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: PetscInitialize(&argc, &argv, (char *)0, help);
76: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: Initialize problem parameters
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: user.param = 6.0;
80: PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create nonlinear solver context
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: SNESCreate(PETSC_COMM_WORLD, &snes);
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Create distributed array (DMDA) to manage parallel grid and vectors
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 4, 4, 4, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &user.da);
92: DMSetFromOptions(user.da);
93: DMSetUp(user.da);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Extract global vectors from DMDA; then duplicate for remaining
96: vectors that are the same types
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: DMCreateGlobalVector(user.da, &x);
99: VecDuplicate(x, &r);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Set function evaluation routine and vector
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: SNESSetFunction(snes, r, FormFunction, (void *)&user);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create matrix data structure; set Jacobian evaluation routine
109: Set Jacobian matrix data structure and default Jacobian evaluation
110: routine. User can override with:
111: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
112: (unless user explicitly sets preconditioner)
113: -snes_mf_operator : form preconditioning matrix as set by the user,
114: but use matrix-free approx for Jacobian-vector
115: products within Newton-Krylov method
116: -fdcoloring : using finite differences with coloring to compute the Jacobian
118: Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option
119: below is to test the call to MatFDColoringSetType().
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: PetscOptionsGetBool(NULL, NULL, "-snes_mf", &matrix_free, NULL);
122: PetscOptionsGetBool(NULL, NULL, "-fdcoloring", &coloring, NULL);
123: PetscOptionsGetBool(NULL, NULL, "-fdcoloring_ds", &coloring_ds, NULL);
124: PetscOptionsGetBool(NULL, NULL, "-fdcoloring_local", &local_coloring, NULL);
125: if (!matrix_free) {
126: DMSetMatType(user.da, MATAIJ);
127: DMCreateMatrix(user.da, &J);
128: if (coloring) {
129: ISColoring iscoloring;
130: if (!local_coloring) {
131: DMCreateColoring(user.da, IS_COLORING_GLOBAL, &iscoloring);
132: MatFDColoringCreate(J, iscoloring, &matfdcoloring);
133: MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunction, &user);
134: } else {
135: DMCreateColoring(user.da, IS_COLORING_LOCAL, &iscoloring);
136: MatFDColoringCreate(J, iscoloring, &matfdcoloring);
137: MatFDColoringUseDM(J, matfdcoloring);
138: MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunctionLocal, &user);
139: }
140: if (coloring_ds) MatFDColoringSetType(matfdcoloring, MATMFFD_DS);
141: MatFDColoringSetFromOptions(matfdcoloring);
142: MatFDColoringSetUp(J, iscoloring, matfdcoloring);
143: SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring);
144: ISColoringDestroy(&iscoloring);
145: } else {
146: SNESSetJacobian(snes, J, J, FormJacobian, &user);
147: }
148: }
150: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: Customize nonlinear solver; set runtime options
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: SNESSetDM(snes, user.da);
154: SNESSetFromOptions(snes);
156: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157: Evaluate initial guess
158: Note: The user should initialize the vector, x, with the initial guess
159: for the nonlinear solver prior to calling SNESSolve(). In particular,
160: to employ an initial guess of zero, the user should explicitly set
161: this vector to zero by calling VecSet().
162: */
163: FormInitialGuess(&user, x);
165: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166: Solve nonlinear system
167: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168: SNESSolve(snes, NULL, x);
169: SNESGetIterationNumber(snes, &its);
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Explicitly check norm of the residual of the solution
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: FormFunction(snes, x, r, (void *)&user);
175: VecNorm(r, NORM_2, &fnorm);
176: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT " fnorm %g\n", its, (double)fnorm);
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Free work space. All PETSc objects should be destroyed when they
180: are no longer needed.
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: MatDestroy(&J);
184: VecDestroy(&x);
185: VecDestroy(&r);
186: SNESDestroy(&snes);
187: DMDestroy(&user.da);
188: MatFDColoringDestroy(&matfdcoloring);
189: PetscFinalize();
190: return 0;
191: }
192: /* ------------------------------------------------------------------- */
193: /*
194: FormInitialGuess - Forms initial approximation.
196: Input Parameters:
197: user - user-defined application context
198: X - vector
200: Output Parameter:
201: X - vector
202: */
203: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
204: {
205: PetscInt i, j, k, Mx, My, Mz, xs, ys, zs, xm, ym, zm;
206: PetscReal lambda, temp1, hx, hy, hz, tempk, tempj;
207: PetscScalar ***x;
210: DMDAGetInfo(user->da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
212: lambda = user->param;
213: hx = 1.0 / (PetscReal)(Mx - 1);
214: hy = 1.0 / (PetscReal)(My - 1);
215: hz = 1.0 / (PetscReal)(Mz - 1);
216: temp1 = lambda / (lambda + 1.0);
218: /*
219: Get a pointer to vector data.
220: - For default PETSc vectors, VecGetArray() returns a pointer to
221: the data array. Otherwise, the routine is implementation dependent.
222: - You MUST call VecRestoreArray() when you no longer need access to
223: the array.
224: */
225: DMDAVecGetArray(user->da, X, &x);
227: /*
228: Get local grid boundaries (for 3-dimensional DMDA):
229: xs, ys, zs - starting grid indices (no ghost points)
230: xm, ym, zm - widths of local grid (no ghost points)
232: */
233: DMDAGetCorners(user->da, &xs, &ys, &zs, &xm, &ym, &zm);
235: /*
236: Compute initial guess over the locally owned part of the grid
237: */
238: for (k = zs; k < zs + zm; k++) {
239: tempk = (PetscReal)(PetscMin(k, Mz - k - 1)) * hz;
240: for (j = ys; j < ys + ym; j++) {
241: tempj = PetscMin((PetscReal)(PetscMin(j, My - j - 1)) * hy, tempk);
242: for (i = xs; i < xs + xm; i++) {
243: if (i == 0 || j == 0 || k == 0 || i == Mx - 1 || j == My - 1 || k == Mz - 1) {
244: /* boundary conditions are all zero Dirichlet */
245: x[k][j][i] = 0.0;
246: } else {
247: x[k][j][i] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, Mx - i - 1)) * hx, tempj));
248: }
249: }
250: }
251: }
253: /*
254: Restore vector
255: */
256: DMDAVecRestoreArray(user->da, X, &x);
257: return 0;
258: }
259: /* ------------------------------------------------------------------- */
260: /*
261: FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch
263: Input Parameters:
264: . snes - the SNES context
265: . localX - input vector, this contains the ghosted region needed
266: . ptr - optional user-defined context, as set by SNESSetFunction()
268: Output Parameter:
269: . F - function vector, this does not contain a ghosted region
270: */
271: PetscErrorCode FormFunctionLocal(SNES snes, Vec localX, Vec F, void *ptr)
272: {
273: AppCtx *user = (AppCtx *)ptr;
274: PetscInt i, j, k, Mx, My, Mz, xs, ys, zs, xm, ym, zm;
275: PetscReal two = 2.0, lambda, hx, hy, hz, hxhzdhy, hyhzdhx, hxhydhz, sc;
276: PetscScalar u_north, u_south, u_east, u_west, u_up, u_down, u, u_xx, u_yy, u_zz, ***x, ***f;
277: DM da;
280: SNESGetDM(snes, &da);
281: DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
283: lambda = user->param;
284: hx = 1.0 / (PetscReal)(Mx - 1);
285: hy = 1.0 / (PetscReal)(My - 1);
286: hz = 1.0 / (PetscReal)(Mz - 1);
287: sc = hx * hy * hz * lambda;
288: hxhzdhy = hx * hz / hy;
289: hyhzdhx = hy * hz / hx;
290: hxhydhz = hx * hy / hz;
292: /*
293: Get pointers to vector data
294: */
295: DMDAVecGetArrayRead(da, localX, &x);
296: DMDAVecGetArray(da, F, &f);
298: /*
299: Get local grid boundaries
300: */
301: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
303: /*
304: Compute function over the locally owned part of the grid
305: */
306: for (k = zs; k < zs + zm; k++) {
307: for (j = ys; j < ys + ym; j++) {
308: for (i = xs; i < xs + xm; i++) {
309: if (i == 0 || j == 0 || k == 0 || i == Mx - 1 || j == My - 1 || k == Mz - 1) {
310: f[k][j][i] = x[k][j][i];
311: } else {
312: u = x[k][j][i];
313: u_east = x[k][j][i + 1];
314: u_west = x[k][j][i - 1];
315: u_north = x[k][j + 1][i];
316: u_south = x[k][j - 1][i];
317: u_up = x[k + 1][j][i];
318: u_down = x[k - 1][j][i];
319: u_xx = (-u_east + two * u - u_west) * hyhzdhx;
320: u_yy = (-u_north + two * u - u_south) * hxhzdhy;
321: u_zz = (-u_up + two * u - u_down) * hxhydhz;
322: f[k][j][i] = u_xx + u_yy + u_zz - sc * PetscExpScalar(u);
323: }
324: }
325: }
326: }
328: /*
329: Restore vectors
330: */
331: DMDAVecRestoreArrayRead(da, localX, &x);
332: DMDAVecRestoreArray(da, F, &f);
333: PetscLogFlops(11.0 * ym * xm);
334: return 0;
335: }
336: /* ------------------------------------------------------------------- */
337: /*
338: FormFunction - Evaluates nonlinear function, F(x) on the entire domain
340: Input Parameters:
341: . snes - the SNES context
342: . X - input vector
343: . ptr - optional user-defined context, as set by SNESSetFunction()
345: Output Parameter:
346: . F - function vector
347: */
348: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
349: {
350: Vec localX;
351: DM da;
354: SNESGetDM(snes, &da);
355: DMGetLocalVector(da, &localX);
357: /*
358: Scatter ghost points to local vector,using the 2-step process
359: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
360: By placing code between these two statements, computations can be
361: done while messages are in transition.
362: */
363: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
364: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
366: FormFunctionLocal(snes, localX, F, ptr);
367: DMRestoreLocalVector(da, &localX);
368: return 0;
369: }
370: /* ------------------------------------------------------------------- */
371: /*
372: FormJacobian - Evaluates Jacobian matrix.
374: Input Parameters:
375: . snes - the SNES context
376: . x - input vector
377: . ptr - optional user-defined context, as set by SNESSetJacobian()
379: Output Parameters:
380: . A - Jacobian matrix
381: . B - optionally different preconditioning matrix
383: */
384: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
385: {
386: AppCtx *user = (AppCtx *)ptr; /* user-defined application context */
387: Vec localX;
388: PetscInt i, j, k, Mx, My, Mz;
389: MatStencil col[7], row;
390: PetscInt xs, ys, zs, xm, ym, zm;
391: PetscScalar lambda, v[7], hx, hy, hz, hxhzdhy, hyhzdhx, hxhydhz, sc, ***x;
392: DM da;
395: SNESGetDM(snes, &da);
396: DMGetLocalVector(da, &localX);
397: DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, &Mz, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE);
399: lambda = user->param;
400: hx = 1.0 / (PetscReal)(Mx - 1);
401: hy = 1.0 / (PetscReal)(My - 1);
402: hz = 1.0 / (PetscReal)(Mz - 1);
403: sc = hx * hy * hz * lambda;
404: hxhzdhy = hx * hz / hy;
405: hyhzdhx = hy * hz / hx;
406: hxhydhz = hx * hy / hz;
408: /*
409: Scatter ghost points to local vector, using the 2-step process
410: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
411: By placing code between these two statements, computations can be
412: done while messages are in transition.
413: */
414: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
415: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
417: /*
418: Get pointer to vector data
419: */
420: DMDAVecGetArrayRead(da, localX, &x);
422: /*
423: Get local grid boundaries
424: */
425: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
427: /*
428: Compute entries for the locally owned part of the Jacobian.
429: - Currently, all PETSc parallel matrix formats are partitioned by
430: contiguous chunks of rows across the processors.
431: - Each processor needs to insert only elements that it owns
432: locally (but any non-local elements will be sent to the
433: appropriate processor during matrix assembly).
434: - Here, we set all entries for a particular row at once.
435: - We can set matrix entries either using either
436: MatSetValuesLocal() or MatSetValues(), as discussed above.
437: */
438: for (k = zs; k < zs + zm; k++) {
439: for (j = ys; j < ys + ym; j++) {
440: for (i = xs; i < xs + xm; i++) {
441: row.k = k;
442: row.j = j;
443: row.i = i;
444: /* boundary points */
445: if (i == 0 || j == 0 || k == 0 || i == Mx - 1 || j == My - 1 || k == Mz - 1) {
446: v[0] = 1.0;
447: MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES);
448: } else {
449: /* interior grid points */
450: v[0] = -hxhydhz;
451: col[0].k = k - 1;
452: col[0].j = j;
453: col[0].i = i;
454: v[1] = -hxhzdhy;
455: col[1].k = k;
456: col[1].j = j - 1;
457: col[1].i = i;
458: v[2] = -hyhzdhx;
459: col[2].k = k;
460: col[2].j = j;
461: col[2].i = i - 1;
462: v[3] = 2.0 * (hyhzdhx + hxhzdhy + hxhydhz) - sc * PetscExpScalar(x[k][j][i]);
463: col[3].k = row.k;
464: col[3].j = row.j;
465: col[3].i = row.i;
466: v[4] = -hyhzdhx;
467: col[4].k = k;
468: col[4].j = j;
469: col[4].i = i + 1;
470: v[5] = -hxhzdhy;
471: col[5].k = k;
472: col[5].j = j + 1;
473: col[5].i = i;
474: v[6] = -hxhydhz;
475: col[6].k = k + 1;
476: col[6].j = j;
477: col[6].i = i;
478: MatSetValuesStencil(jac, 1, &row, 7, col, v, INSERT_VALUES);
479: }
480: }
481: }
482: }
483: DMDAVecRestoreArrayRead(da, localX, &x);
484: DMRestoreLocalVector(da, &localX);
486: /*
487: Assemble matrix, using the 2-step process:
488: MatAssemblyBegin(), MatAssemblyEnd().
489: */
490: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
491: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
493: /*
494: Normally since the matrix has already been assembled above; this
495: would do nothing. But in the matrix free mode -snes_mf_operator
496: this tells the "matrix-free" matrix that a new linear system solve
497: is about to be done.
498: */
500: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
501: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
503: /*
504: Tell the matrix we will never add a new nonzero location to the
505: matrix. If we do, it will generate an error.
506: */
507: MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
508: return 0;
509: }
511: /*TEST
513: test:
514: nsize: 4
515: args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
517: test:
518: suffix: 2
519: nsize: 4
520: args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
522: test:
523: suffix: 3
524: nsize: 4
525: args: -fdcoloring -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
527: test:
528: suffix: 3_ds
529: nsize: 4
530: args: -fdcoloring -fdcoloring_ds -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
532: test:
533: suffix: 4
534: nsize: 4
535: args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1
536: requires: !single
538: test:
539: suffix: 5
540: nsize: 4
541: args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1 -snes_type newtontrdc
542: requires: !single
544: TEST*/