Actual source code: ex14.c
2: static char help[] = "Solves a nonlinear system in parallel with a user-defined Newton method.\n\
3: Uses KSP to solve the linearized Newton systems. This solver\n\
4: is a very simplistic inexact Newton method. The intent of this code is to\n\
5: demonstrate the repeated solution of linear systems with the same nonzero pattern.\n\
6: \n\
7: This is NOT the recommended approach for solving nonlinear problems with PETSc!\n\
8: We urge users to employ the SNES component for solving nonlinear problems whenever\n\
9: possible, as it offers many advantages over coding nonlinear solvers independently.\n\
10: \n\
11: We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
12: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
13: The command line options include:\n\
14: -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
15: problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
16: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
17: -my <yg>, where <yg> = number of grid points in the y-direction\n\
18: -Nx <npx>, where <npx> = number of processors in the x-direction\n\
19: -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";
21: /* ------------------------------------------------------------------------
23: Solid Fuel Ignition (SFI) problem. This problem is modeled by
24: the partial differential equation
26: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
28: with boundary conditions
30: u = 0 for x = 0, x = 1, y = 0, y = 1.
32: A finite difference approximation with the usual 5-point stencil
33: is used to discretize the boundary value problem to obtain a nonlinear
34: system of equations.
36: The SNES version of this problem is: snes/tutorials/ex5.c
37: We urge users to employ the SNES component for solving nonlinear
38: problems whenever possible, as it offers many advantages over coding
39: nonlinear solvers independently.
41: ------------------------------------------------------------------------- */
43: /*
44: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
45: Include "petscksp.h" so that we can use KSP solvers. Note that this
46: file automatically includes:
47: petscsys.h - base PETSc routines petscvec.h - vectors
48: petscmat.h - matrices
49: petscis.h - index sets petscksp.h - Krylov subspace methods
50: petscviewer.h - viewers petscpc.h - preconditioners
51: */
52: #include <petscdm.h>
53: #include <petscdmda.h>
54: #include <petscksp.h>
56: /*
57: User-defined application context - contains data needed by the
58: application-provided call-back routines, ComputeJacobian() and
59: ComputeFunction().
60: */
61: typedef struct {
62: PetscReal param; /* test problem parameter */
63: PetscInt mx, my; /* discretization in x,y directions */
64: Vec localX; /* ghosted local vector */
65: DM da; /* distributed array data structure */
66: } AppCtx;
68: /*
69: User-defined routines
70: */
71: extern PetscErrorCode ComputeFunction(AppCtx *, Vec, Vec), FormInitialGuess(AppCtx *, Vec);
72: extern PetscErrorCode ComputeJacobian(AppCtx *, Vec, Mat);
74: int main(int argc, char **argv)
75: {
76: /* -------------- Data to define application problem ---------------- */
77: MPI_Comm comm; /* communicator */
78: KSP ksp; /* linear solver */
79: Vec X, Y, F; /* solution, update, residual vectors */
80: Mat J; /* Jacobian matrix */
81: AppCtx user; /* user-defined work context */
82: PetscInt Nx, Ny; /* number of preocessors in x- and y- directions */
83: PetscMPIInt size; /* number of processors */
84: PetscReal bratu_lambda_max = 6.81, bratu_lambda_min = 0.;
85: PetscInt m, N;
87: /* --------------- Data to define nonlinear solver -------------- */
88: PetscReal rtol = 1.e-8; /* relative convergence tolerance */
89: PetscReal xtol = 1.e-8; /* step convergence tolerance */
90: PetscReal ttol; /* convergence tolerance */
91: PetscReal fnorm, ynorm, xnorm; /* various vector norms */
92: PetscInt max_nonlin_its = 3; /* maximum number of iterations for nonlinear solver */
93: PetscInt max_functions = 50; /* maximum number of function evaluations */
94: PetscInt lin_its; /* number of linear solver iterations for each step */
95: PetscInt i; /* nonlinear solve iteration number */
96: PetscBool no_output = PETSC_FALSE; /* flag indicating whether to suppress output */
99: PetscInitialize(&argc, &argv, (char *)0, help);
100: comm = PETSC_COMM_WORLD;
101: PetscOptionsGetBool(NULL, NULL, "-no_output", &no_output, NULL);
103: /*
104: Initialize problem parameters
105: */
106: user.mx = 4;
107: user.my = 4;
108: user.param = 6.0;
110: PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL);
111: PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL);
112: PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL);
114: N = user.mx * user.my;
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Create linear solver context
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: KSPCreate(comm, &ksp);
122: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: Create vector data structures
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: /*
127: Create distributed array (DMDA) to manage parallel grid and vectors
128: */
129: MPI_Comm_size(comm, &size);
130: Nx = PETSC_DECIDE;
131: Ny = PETSC_DECIDE;
132: PetscOptionsGetInt(NULL, NULL, "-Nx", &Nx, NULL);
133: PetscOptionsGetInt(NULL, NULL, "-Ny", &Ny, NULL);
135: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, user.mx, user.my, Nx, Ny, 1, 1, NULL, NULL, &user.da);
136: DMSetFromOptions(user.da);
137: DMSetUp(user.da);
139: /*
140: Extract global and local vectors from DMDA; then duplicate for remaining
141: vectors that are the same types
142: */
143: DMCreateGlobalVector(user.da, &X);
144: DMCreateLocalVector(user.da, &user.localX);
145: VecDuplicate(X, &F);
146: VecDuplicate(X, &Y);
148: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: Create matrix data structure for Jacobian
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151: /*
152: Note: For the parallel case, vectors and matrices MUST be partitioned
153: accordingly. When using distributed arrays (DMDAs) to create vectors,
154: the DMDAs determine the problem partitioning. We must explicitly
155: specify the local matrix dimensions upon its creation for compatibility
156: with the vector distribution. Thus, the generic MatCreate() routine
157: is NOT sufficient when working with distributed arrays.
159: Note: Here we only approximately preallocate storage space for the
160: Jacobian. See the users manual for a discussion of better techniques
161: for preallocating matrix memory.
162: */
163: if (size == 1) {
164: MatCreateSeqAIJ(comm, N, N, 5, NULL, &J);
165: } else {
166: VecGetLocalSize(X, &m);
167: MatCreateAIJ(comm, m, m, N, N, 5, NULL, 3, NULL, &J);
168: }
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Customize linear solver; set runtime options
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: /*
175: Set runtime options (e.g.,-ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
176: */
177: KSPSetFromOptions(ksp);
179: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: Evaluate initial guess
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: FormInitialGuess(&user, X);
184: ComputeFunction(&user, X, F); /* Compute F(X) */
185: VecNorm(F, NORM_2, &fnorm); /* fnorm = || F || */
186: ttol = fnorm * rtol;
187: if (!no_output) PetscPrintf(comm, "Initial function norm = %g\n", (double)fnorm);
189: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: Solve nonlinear system with a user-defined method
191: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
193: /*
194: This solver is a very simplistic inexact Newton method, with no
195: no damping strategies or bells and whistles. The intent of this code
196: is merely to demonstrate the repeated solution with KSP of linear
197: systems with the same nonzero structure.
199: This is NOT the recommended approach for solving nonlinear problems
200: with PETSc! We urge users to employ the SNES component for solving
201: nonlinear problems whenever possible with application codes, as it
202: offers many advantages over coding nonlinear solvers independently.
203: */
205: for (i = 0; i < max_nonlin_its; i++) {
206: /*
207: Compute the Jacobian matrix.
208: */
209: ComputeJacobian(&user, X, J);
211: /*
212: Solve J Y = F, where J is the Jacobian matrix.
213: - First, set the KSP linear operators. Here the matrix that
214: defines the linear system also serves as the preconditioning
215: matrix.
216: - Then solve the Newton system.
217: */
218: KSPSetOperators(ksp, J, J);
219: KSPSolve(ksp, F, Y);
220: KSPGetIterationNumber(ksp, &lin_its);
222: /*
223: Compute updated iterate
224: */
225: VecNorm(Y, NORM_2, &ynorm); /* ynorm = || Y || */
226: VecAYPX(Y, -1.0, X); /* Y <- X - Y */
227: VecCopy(Y, X); /* X <- Y */
228: VecNorm(X, NORM_2, &xnorm); /* xnorm = || X || */
229: if (!no_output) PetscPrintf(comm, " linear solve iterations = %" PetscInt_FMT ", xnorm=%g, ynorm=%g\n", lin_its, (double)xnorm, (double)ynorm);
231: /*
232: Evaluate new nonlinear function
233: */
234: ComputeFunction(&user, X, F); /* Compute F(X) */
235: VecNorm(F, NORM_2, &fnorm); /* fnorm = || F || */
236: if (!no_output) PetscPrintf(comm, "Iteration %" PetscInt_FMT ", function norm = %g\n", i + 1, (double)fnorm);
238: /*
239: Test for convergence
240: */
241: if (fnorm <= ttol) {
242: if (!no_output) PetscPrintf(comm, "Converged due to function norm %g < %g (relative tolerance)\n", (double)fnorm, (double)ttol);
243: break;
244: }
245: if (ynorm < xtol * (xnorm)) {
246: if (!no_output) PetscPrintf(comm, "Converged due to small update length: %g < %g * %g\n", (double)ynorm, (double)xtol, (double)xnorm);
247: break;
248: }
249: if (i > max_functions) {
250: if (!no_output) PetscPrintf(comm, "Exceeded maximum number of function evaluations: %" PetscInt_FMT " > %" PetscInt_FMT "\n", i, max_functions);
251: break;
252: }
253: }
254: PetscPrintf(comm, "Number of nonlinear iterations = %" PetscInt_FMT "\n", i);
256: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257: Free work space. All PETSc objects should be destroyed when they
258: are no longer needed.
259: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
261: MatDestroy(&J);
262: VecDestroy(&Y);
263: VecDestroy(&user.localX);
264: VecDestroy(&X);
265: VecDestroy(&F);
266: KSPDestroy(&ksp);
267: DMDestroy(&user.da);
268: PetscFinalize();
269: return 0;
270: }
271: /* ------------------------------------------------------------------- */
272: /*
273: FormInitialGuess - Forms initial approximation.
275: Input Parameters:
276: user - user-defined application context
277: X - vector
279: Output Parameter:
280: X - vector
281: */
282: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
283: {
284: PetscInt i, j, row, mx, my, xs, ys, xm, ym, gxm, gym, gxs, gys;
285: PetscReal one = 1.0, lambda, temp1, temp, hx, hy;
286: PetscScalar *x;
288: mx = user->mx;
289: my = user->my;
290: lambda = user->param;
291: hx = one / (PetscReal)(mx - 1);
292: hy = one / (PetscReal)(my - 1);
293: temp1 = lambda / (lambda + one);
295: /*
296: Get a pointer to vector data.
297: - For default PETSc vectors, VecGetArray() returns a pointer to
298: the data array. Otherwise, the routine is implementation dependent.
299: - You MUST call VecRestoreArray() when you no longer need access to
300: the array.
301: */
302: VecGetArray(X, &x);
304: /*
305: Get local grid boundaries (for 2-dimensional DMDA):
306: xs, ys - starting grid indices (no ghost points)
307: xm, ym - widths of local grid (no ghost points)
308: gxs, gys - starting grid indices (including ghost points)
309: gxm, gym - widths of local grid (including ghost points)
310: */
311: DMDAGetCorners(user->da, &xs, &ys, NULL, &xm, &ym, NULL);
312: DMDAGetGhostCorners(user->da, &gxs, &gys, NULL, &gxm, &gym, NULL);
314: /*
315: Compute initial guess over the locally owned part of the grid
316: */
317: for (j = ys; j < ys + ym; j++) {
318: temp = (PetscReal)(PetscMin(j, my - j - 1)) * hy;
319: for (i = xs; i < xs + xm; i++) {
320: row = i - gxs + (j - gys) * gxm;
321: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
322: x[row] = 0.0;
323: continue;
324: }
325: x[row] = temp1 * PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i, mx - i - 1)) * hx, temp));
326: }
327: }
329: /*
330: Restore vector
331: */
332: VecRestoreArray(X, &x);
333: return 0;
334: }
335: /* ------------------------------------------------------------------- */
336: /*
337: ComputeFunction - Evaluates nonlinear function, F(x).
339: Input Parameters:
340: . X - input vector
341: . user - user-defined application context
343: Output Parameter:
344: . F - function vector
345: */
346: PetscErrorCode ComputeFunction(AppCtx *user, Vec X, Vec F)
347: {
348: PetscInt i, j, row, mx, my, xs, ys, xm, ym, gxs, gys, gxm, gym;
349: PetscReal two = 2.0, one = 1.0, lambda, hx, hy, hxdhy, hydhx, sc;
350: PetscScalar u, uxx, uyy, *x, *f;
351: Vec localX = user->localX;
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 * lambda;
359: hxdhy = hx / hy;
360: hydhx = hy / hx;
362: /*
363: Scatter ghost points to local vector, using the 2-step process
364: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
365: By placing code between these two statements, computations can be
366: done while messages are in transition.
367: */
368: DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX);
369: DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX);
371: /*
372: Get pointers to vector data
373: */
374: VecGetArray(localX, &x);
375: VecGetArray(F, &f);
377: /*
378: Get local grid boundaries
379: */
380: DMDAGetCorners(user->da, &xs, &ys, NULL, &xm, &ym, NULL);
381: DMDAGetGhostCorners(user->da, &gxs, &gys, NULL, &gxm, &gym, NULL);
383: /*
384: Compute function over the locally owned part of the grid
385: */
386: for (j = ys; j < ys + ym; j++) {
387: row = (j - gys) * gxm + xs - gxs - 1;
388: for (i = xs; i < xs + xm; i++) {
389: row++;
390: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
391: f[row] = x[row];
392: continue;
393: }
394: u = x[row];
395: uxx = (two * u - x[row - 1] - x[row + 1]) * hydhx;
396: uyy = (two * u - x[row - gxm] - x[row + gxm]) * hxdhy;
397: f[row] = uxx + uyy - sc * PetscExpScalar(u);
398: }
399: }
401: /*
402: Restore vectors
403: */
404: VecRestoreArray(localX, &x);
405: VecRestoreArray(F, &f);
406: PetscLogFlops(11.0 * ym * xm);
407: return 0;
408: }
409: /* ------------------------------------------------------------------- */
410: /*
411: ComputeJacobian - Evaluates Jacobian matrix.
413: Input Parameters:
414: . x - input vector
415: . user - user-defined application context
417: Output Parameters:
418: . jac - Jacobian matrix
419: . flag - flag indicating matrix structure
421: Notes:
422: Due to grid point reordering with DMDAs, we must always work
423: with the local grid points, and then transform them to the new
424: global numbering with the "ltog" mapping
425: We cannot work directly with the global numbers for the original
426: uniprocessor grid!
427: */
428: PetscErrorCode ComputeJacobian(AppCtx *user, Vec X, Mat jac)
429: {
430: Vec localX = user->localX; /* local vector */
431: const PetscInt *ltog; /* local-to-global mapping */
432: PetscInt i, j, row, mx, my, col[5];
433: PetscInt xs, ys, xm, ym, gxs, gys, gxm, gym, grow;
434: PetscScalar two = 2.0, one = 1.0, lambda, v[5], hx, hy, hxdhy, hydhx, sc, *x;
435: ISLocalToGlobalMapping ltogm;
437: mx = user->mx;
438: my = user->my;
439: lambda = user->param;
440: hx = one / (PetscReal)(mx - 1);
441: hy = one / (PetscReal)(my - 1);
442: sc = hx * hy;
443: hxdhy = hx / hy;
444: hydhx = hy / hx;
446: /*
447: Scatter ghost points to local vector, using the 2-step process
448: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
449: By placing code between these two statements, computations can be
450: done while messages are in transition.
451: */
452: DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX);
453: DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX);
455: /*
456: Get pointer to vector data
457: */
458: VecGetArray(localX, &x);
460: /*
461: Get local grid boundaries
462: */
463: DMDAGetCorners(user->da, &xs, &ys, NULL, &xm, &ym, NULL);
464: DMDAGetGhostCorners(user->da, &gxs, &gys, NULL, &gxm, &gym, NULL);
466: /*
467: Get the global node numbers for all local nodes, including ghost points
468: */
469: DMGetLocalToGlobalMapping(user->da, <ogm);
470: ISLocalToGlobalMappingGetIndices(ltogm, <og);
472: /*
473: Compute entries for the locally owned part of the Jacobian.
474: - Currently, all PETSc parallel matrix formats are partitioned by
475: contiguous chunks of rows across the processors. The "grow"
476: parameter computed below specifies the global row number
477: corresponding to each local grid point.
478: - Each processor needs to insert only elements that it owns
479: locally (but any non-local elements will be sent to the
480: appropriate processor during matrix assembly).
481: - Always specify global row and columns of matrix entries.
482: - Here, we set all entries for a particular row at once.
483: */
484: for (j = ys; j < ys + ym; j++) {
485: row = (j - gys) * gxm + xs - gxs - 1;
486: for (i = xs; i < xs + xm; i++) {
487: row++;
488: grow = ltog[row];
489: /* boundary points */
490: if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
491: MatSetValues(jac, 1, &grow, 1, &grow, &one, INSERT_VALUES);
492: continue;
493: }
494: /* interior grid points */
495: v[0] = -hxdhy;
496: col[0] = ltog[row - gxm];
497: v[1] = -hydhx;
498: col[1] = ltog[row - 1];
499: v[2] = two * (hydhx + hxdhy) - sc * lambda * PetscExpScalar(x[row]);
500: col[2] = grow;
501: v[3] = -hydhx;
502: col[3] = ltog[row + 1];
503: v[4] = -hxdhy;
504: col[4] = ltog[row + gxm];
505: MatSetValues(jac, 1, &grow, 5, col, v, INSERT_VALUES);
506: }
507: }
508: ISLocalToGlobalMappingRestoreIndices(ltogm, <og);
510: /*
511: Assemble matrix, using the 2-step process:
512: MatAssemblyBegin(), MatAssemblyEnd().
513: By placing code between these two statements, computations can be
514: done while messages are in transition.
515: */
516: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
517: VecRestoreArray(localX, &x);
518: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
520: return 0;
521: }
523: /*TEST
525: test:
527: TEST*/