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 surpress output */
98: PetscInitialize(&argc,&argv,(char*)0,help);
99: comm = PETSC_COMM_WORLD;
100: PetscOptionsGetBool(NULL,NULL,"-no_output",&no_output,NULL);
102: /*
103: Initialize problem parameters
104: */
105: user.mx = 4; user.my = 4; user.param = 6.0;
107: PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
108: PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
109: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
111: N = user.mx*user.my;
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Create linear solver context
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: KSPCreate(comm,&ksp);
119: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120: Create vector data structures
121: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123: /*
124: Create distributed array (DMDA) to manage parallel grid and vectors
125: */
126: MPI_Comm_size(comm,&size);
127: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
128: PetscOptionsGetInt(NULL,NULL,"-Nx",&Nx,NULL);
129: PetscOptionsGetInt(NULL,NULL,"-Ny",&Ny,NULL);
131: 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);
132: DMSetFromOptions(user.da);
133: DMSetUp(user.da);
135: /*
136: Extract global and local vectors from DMDA; then duplicate for remaining
137: vectors that are the same types
138: */
139: DMCreateGlobalVector(user.da,&X);
140: DMCreateLocalVector(user.da,&user.localX);
141: VecDuplicate(X,&F);
142: VecDuplicate(X,&Y);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Create matrix data structure for Jacobian
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: /*
148: Note: For the parallel case, vectors and matrices MUST be partitioned
149: accordingly. When using distributed arrays (DMDAs) to create vectors,
150: the DMDAs determine the problem partitioning. We must explicitly
151: specify the local matrix dimensions upon its creation for compatibility
152: with the vector distribution. Thus, the generic MatCreate() routine
153: is NOT sufficient when working with distributed arrays.
155: Note: Here we only approximately preallocate storage space for the
156: Jacobian. See the users manual for a discussion of better techniques
157: for preallocating matrix memory.
158: */
159: if (size == 1) {
160: MatCreateSeqAIJ(comm,N,N,5,NULL,&J);
161: } else {
162: VecGetLocalSize(X,&m);
163: MatCreateAIJ(comm,m,m,N,N,5,NULL,3,NULL,&J);
164: }
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Customize linear solver; set runtime options
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170: /*
171: Set runtime options (e.g.,-ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
172: */
173: KSPSetFromOptions(ksp);
175: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176: Evaluate initial guess
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: FormInitialGuess(&user,X);
180: ComputeFunction(&user,X,F); /* Compute F(X) */
181: VecNorm(F,NORM_2,&fnorm); /* fnorm = || F || */
182: ttol = fnorm*rtol;
183: if (!no_output) PetscPrintf(comm,"Initial function norm = %g\n",(double)fnorm);
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Solve nonlinear system with a user-defined method
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: /*
190: This solver is a very simplistic inexact Newton method, with no
191: no damping strategies or bells and whistles. The intent of this code
192: is merely to demonstrate the repeated solution with KSP of linear
193: systems with the same nonzero structure.
195: This is NOT the recommended approach for solving nonlinear problems
196: with PETSc! We urge users to employ the SNES component for solving
197: nonlinear problems whenever possible with application codes, as it
198: offers many advantages over coding nonlinear solvers independently.
199: */
201: for (i=0; i<max_nonlin_its; i++) {
203: /*
204: Compute the Jacobian matrix.
205: */
206: ComputeJacobian(&user,X,J);
208: /*
209: Solve J Y = F, where J is the Jacobian matrix.
210: - First, set the KSP linear operators. Here the matrix that
211: defines the linear system also serves as the preconditioning
212: matrix.
213: - Then solve the Newton system.
214: */
215: KSPSetOperators(ksp,J,J);
216: KSPSolve(ksp,F,Y);
217: KSPGetIterationNumber(ksp,&lin_its);
219: /*
220: Compute updated iterate
221: */
222: VecNorm(Y,NORM_2,&ynorm); /* ynorm = || Y || */
223: VecAYPX(Y,-1.0,X); /* Y <- X - Y */
224: VecCopy(Y,X); /* X <- Y */
225: VecNorm(X,NORM_2,&xnorm); /* xnorm = || X || */
226: if (!no_output) {
227: PetscPrintf(comm," linear solve iterations = %D, xnorm=%g, ynorm=%g\n",lin_its,(double)xnorm,(double)ynorm);
228: }
230: /*
231: Evaluate new nonlinear function
232: */
233: ComputeFunction(&user,X,F); /* Compute F(X) */
234: VecNorm(F,NORM_2,&fnorm); /* fnorm = || F || */
235: if (!no_output) {
236: PetscPrintf(comm,"Iteration %D, function norm = %g\n",i+1,(double)fnorm);
237: }
239: /*
240: Test for convergence
241: */
242: if (fnorm <= ttol) {
243: if (!no_output) {
244: PetscPrintf(comm,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)ttol);
245: }
246: break;
247: }
248: if (ynorm < xtol*(xnorm)) {
249: if (!no_output) {
250: PetscPrintf(comm,"Converged due to small update length: %g < %g * %g\n",(double)ynorm,(double)xtol,(double)xnorm);
251: }
252: break;
253: }
254: if (i > max_functions) {
255: if (!no_output) {
256: PetscPrintf(comm,"Exceeded maximum number of function evaluations: %D > %D\n",i,max_functions);
257: }
258: break;
259: }
260: }
261: PetscPrintf(comm,"Number of nonlinear iterations = %D\n",i);
263: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264: Free work space. All PETSc objects should be destroyed when they
265: are no longer needed.
266: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
268: MatDestroy(&J)); PetscCall(VecDestroy(&Y);
269: VecDestroy(&user.localX)); PetscCall(VecDestroy(&X);
270: VecDestroy(&F);
271: KSPDestroy(&ksp)); PetscCall(DMDestroy(&user.da);
272: PetscFinalize();
273: return 0;
274: }
275: /* ------------------------------------------------------------------- */
276: /*
277: FormInitialGuess - Forms initial approximation.
279: Input Parameters:
280: user - user-defined application context
281: X - vector
283: Output Parameter:
284: X - vector
285: */
286: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
287: {
288: PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxm,gym,gxs,gys;
289: PetscReal one = 1.0,lambda,temp1,temp,hx,hy;
290: PetscScalar *x;
292: mx = user->mx; my = user->my; lambda = user->param;
293: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
294: temp1 = lambda/(lambda + one);
296: /*
297: Get a pointer to vector data.
298: - For default PETSc vectors, VecGetArray() returns a pointer to
299: the data array. Otherwise, the routine is implementation dependent.
300: - You MUST call VecRestoreArray() when you no longer need access to
301: the array.
302: */
303: VecGetArray(X,&x);
305: /*
306: Get local grid boundaries (for 2-dimensional DMDA):
307: xs, ys - starting grid indices (no ghost points)
308: xm, ym - widths of local grid (no ghost points)
309: gxs, gys - starting grid indices (including ghost points)
310: gxm, gym - widths of local grid (including ghost points)
311: */
312: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
313: DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);
315: /*
316: Compute initial guess over the locally owned part of the grid
317: */
318: for (j=ys; j<ys+ym; j++) {
319: temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
320: for (i=xs; i<xs+xm; i++) {
321: row = i - gxs + (j - gys)*gxm;
322: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
323: x[row] = 0.0;
324: continue;
325: }
326: x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
327: }
328: }
330: /*
331: Restore vector
332: */
333: VecRestoreArray(X,&x);
334: return 0;
335: }
336: /* ------------------------------------------------------------------- */
337: /*
338: ComputeFunction - Evaluates nonlinear function, F(x).
340: Input Parameters:
341: . X - input vector
342: . user - user-defined application context
344: Output Parameter:
345: . F - function vector
346: */
347: PetscErrorCode ComputeFunction(AppCtx *user,Vec X,Vec F)
348: {
349: PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
350: PetscReal two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
351: PetscScalar u,uxx,uyy,*x,*f;
352: Vec localX = user->localX;
354: mx = user->mx; my = user->my; lambda = user->param;
355: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
356: sc = hx*hy*lambda; hxdhy = hx/hy; hydhx = hy/hx;
358: /*
359: Scatter ghost points to local vector, using the 2-step process
360: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
361: By placing code between these two statements, computations can be
362: done while messages are in transition.
363: */
364: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
365: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
367: /*
368: Get pointers to vector data
369: */
370: VecGetArray(localX,&x);
371: VecGetArray(F,&f);
373: /*
374: Get local grid boundaries
375: */
376: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
377: DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);
379: /*
380: Compute function over the locally owned part of the grid
381: */
382: for (j=ys; j<ys+ym; j++) {
383: row = (j - gys)*gxm + xs - gxs - 1;
384: for (i=xs; i<xs+xm; i++) {
385: row++;
386: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
387: f[row] = x[row];
388: continue;
389: }
390: u = x[row];
391: uxx = (two*u - x[row-1] - x[row+1])*hydhx;
392: uyy = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
393: f[row] = uxx + uyy - sc*PetscExpScalar(u);
394: }
395: }
397: /*
398: Restore vectors
399: */
400: VecRestoreArray(localX,&x);
401: VecRestoreArray(F,&f);
402: PetscLogFlops(11.0*ym*xm);
403: return 0;
404: }
405: /* ------------------------------------------------------------------- */
406: /*
407: ComputeJacobian - Evaluates Jacobian matrix.
409: Input Parameters:
410: . x - input vector
411: . user - user-defined application context
413: Output Parameters:
414: . jac - Jacobian matrix
415: . flag - flag indicating matrix structure
417: Notes:
418: Due to grid point reordering with DMDAs, we must always work
419: with the local grid points, and then transform them to the new
420: global numbering with the "ltog" mapping
421: We cannot work directly with the global numbers for the original
422: uniprocessor grid!
423: */
424: PetscErrorCode ComputeJacobian(AppCtx *user,Vec X,Mat jac)
425: {
426: Vec localX = user->localX; /* local vector */
427: const PetscInt *ltog; /* local-to-global mapping */
428: PetscInt i,j,row,mx,my,col[5];
429: PetscInt xs,ys,xm,ym,gxs,gys,gxm,gym,grow;
430: PetscScalar two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;
431: ISLocalToGlobalMapping ltogm;
433: mx = user->mx; my = user->my; lambda = user->param;
434: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
435: sc = hx*hy; hxdhy = hx/hy; hydhx = hy/hx;
437: /*
438: Scatter ghost points to local vector, using the 2-step process
439: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
440: By placing code between these two statements, computations can be
441: done while messages are in transition.
442: */
443: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
444: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
446: /*
447: Get pointer to vector data
448: */
449: VecGetArray(localX,&x);
451: /*
452: Get local grid boundaries
453: */
454: DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
455: DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);
457: /*
458: Get the global node numbers for all local nodes, including ghost points
459: */
460: DMGetLocalToGlobalMapping(user->da,<ogm);
461: ISLocalToGlobalMappingGetIndices(ltogm,<og);
463: /*
464: Compute entries for the locally owned part of the Jacobian.
465: - Currently, all PETSc parallel matrix formats are partitioned by
466: contiguous chunks of rows across the processors. The "grow"
467: parameter computed below specifies the global row number
468: corresponding to each local grid point.
469: - Each processor needs to insert only elements that it owns
470: locally (but any non-local elements will be sent to the
471: appropriate processor during matrix assembly).
472: - Always specify global row and columns of matrix entries.
473: - Here, we set all entries for a particular row at once.
474: */
475: for (j=ys; j<ys+ym; j++) {
476: row = (j - gys)*gxm + xs - gxs - 1;
477: for (i=xs; i<xs+xm; i++) {
478: row++;
479: grow = ltog[row];
480: /* boundary points */
481: if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
482: MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
483: continue;
484: }
485: /* interior grid points */
486: v[0] = -hxdhy; col[0] = ltog[row - gxm];
487: v[1] = -hydhx; col[1] = ltog[row - 1];
488: v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
489: v[3] = -hydhx; col[3] = ltog[row + 1];
490: v[4] = -hxdhy; col[4] = ltog[row + gxm];
491: MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
492: }
493: }
494: ISLocalToGlobalMappingRestoreIndices(ltogm,<og);
496: /*
497: Assemble matrix, using the 2-step process:
498: MatAssemblyBegin(), MatAssemblyEnd().
499: By placing code between these two statements, computations can be
500: done while messages are in transition.
501: */
502: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
503: VecRestoreArray(localX,&x);
504: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
506: return 0;
507: }
509: /*TEST
511: test:
513: TEST*/