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: /*T
10: Concepts: SNES^parallel Bratu example
11: Concepts: DMDA^using distributed arrays;
12: Processors: n
13: T*/
15: /* ------------------------------------------------------------------------
17: Solid Fuel Ignition (SFI) problem. This problem is modeled by
18: the partial differential equation
20: -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1,
22: with boundary conditions
24: u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
26: A finite difference approximation with the usual 7-point stencil
27: is used to discretize the boundary value problem to obtain a nonlinear
28: system of equations.
30: ------------------------------------------------------------------------- */
32: /*
33: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
34: Include "petscsnes.h" so that we can use SNES solvers. Note that this
35: file automatically includes:
36: petscsys.h - base PETSc routines petscvec.h - vectors
37: petscmat.h - matrices
38: petscis.h - index sets petscksp.h - Krylov subspace methods
39: petscviewer.h - viewers petscpc.h - preconditioners
40: petscksp.h - linear solvers
41: */
42: #include <petscdm.h>
43: #include <petscdmda.h>
44: #include <petscsnes.h>
46: /*
47: User-defined application context - contains data needed by the
48: application-provided call-back routines, FormJacobian() and
49: FormFunction().
50: */
51: typedef struct {
52: PetscReal param; /* test problem parameter */
53: DM da; /* distributed array data structure */
54: } AppCtx;
56: /*
57: User-defined routines
58: */
59: extern PetscErrorCode FormFunctionLocal(SNES,Vec,Vec,void*);
60: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
61: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
62: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
64: int main(int argc,char **argv)
65: {
66: SNES snes; /* nonlinear solver */
67: Vec x,r; /* solution, residual vectors */
68: Mat J = NULL; /* Jacobian matrix */
69: AppCtx user; /* user-defined work context */
70: PetscInt its; /* iterations for convergence */
71: MatFDColoring matfdcoloring = NULL;
72: PetscBool matrix_free = PETSC_FALSE,coloring = PETSC_FALSE, coloring_ds = PETSC_FALSE,local_coloring = PETSC_FALSE;
74: PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm;
76: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: Initialize program
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Initialize problem parameters
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: user.param = 6.0;
86: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
87: if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lambda is out of range");
89: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
90: Create nonlinear solver context
91: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: SNESCreate(PETSC_COMM_WORLD,&snes);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create distributed array (DMDA) to manage parallel grid and vectors
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: 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);
98: DMSetFromOptions(user.da);
99: DMSetUp(user.da);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Extract global vectors from DMDA; then duplicate for remaining
102: vectors that are the same types
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: DMCreateGlobalVector(user.da,&x);
105: VecDuplicate(x,&r);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Set function evaluation routine and vector
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: SNESSetFunction(snes,r,FormFunction,(void*)&user);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create matrix data structure; set Jacobian evaluation routine
115: Set Jacobian matrix data structure and default Jacobian evaluation
116: routine. User can override with:
117: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
118: (unless user explicitly sets preconditioner)
119: -snes_mf_operator : form preconditioning matrix as set by the user,
120: but use matrix-free approx for Jacobian-vector
121: products within Newton-Krylov method
122: -fdcoloring : using finite differences with coloring to compute the Jacobian
124: Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option
125: below is to test the call to MatFDColoringSetType().
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
128: PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL);
129: PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL);
130: PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL);
131: if (!matrix_free) {
132: DMSetMatType(user.da,MATAIJ);
133: DMCreateMatrix(user.da,&J);
134: if (coloring) {
135: ISColoring iscoloring;
136: if (!local_coloring) {
137: DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);
138: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
139: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
140: } else {
141: DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring);
142: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
143: MatFDColoringUseDM(J,matfdcoloring);
144: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunctionLocal,&user);
145: }
146: if (coloring_ds) {
147: MatFDColoringSetType(matfdcoloring,MATMFFD_DS);
148: }
149: MatFDColoringSetFromOptions(matfdcoloring);
150: MatFDColoringSetUp(J,iscoloring,matfdcoloring);
151: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
152: ISColoringDestroy(&iscoloring);
153: } else {
154: SNESSetJacobian(snes,J,J,FormJacobian,&user);
155: }
156: }
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Customize nonlinear solver; set runtime options
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: SNESSetDM(snes,user.da);
162: SNESSetFromOptions(snes);
164: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165: Evaluate initial guess
166: Note: The user should initialize the vector, x, with the initial guess
167: for the nonlinear solver prior to calling SNESSolve(). In particular,
168: to employ an initial guess of zero, the user should explicitly set
169: this vector to zero by calling VecSet().
170: */
171: FormInitialGuess(&user,x);
173: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174: Solve nonlinear system
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: SNESSolve(snes,NULL,x);
177: SNESGetIterationNumber(snes,&its);
179: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: Explicitly check norm of the residual of the solution
181: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: FormFunction(snes,x,r,(void*)&user);
183: VecNorm(r,NORM_2,&fnorm);
184: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm);
186: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
187: Free work space. All PETSc objects should be destroyed when they
188: are no longer needed.
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: MatDestroy(&J);
192: VecDestroy(&x);
193: VecDestroy(&r);
194: SNESDestroy(&snes);
195: DMDestroy(&user.da);
196: MatFDColoringDestroy(&matfdcoloring);
197: PetscFinalize();
198: return ierr;
199: }
200: /* ------------------------------------------------------------------- */
201: /*
202: FormInitialGuess - Forms initial approximation.
204: Input Parameters:
205: user - user-defined application context
206: X - vector
208: Output Parameter:
209: X - vector
210: */
211: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
212: {
213: PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
215: PetscReal lambda,temp1,hx,hy,hz,tempk,tempj;
216: PetscScalar ***x;
219: 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);
221: lambda = user->param;
222: hx = 1.0/(PetscReal)(Mx-1);
223: hy = 1.0/(PetscReal)(My-1);
224: hz = 1.0/(PetscReal)(Mz-1);
225: temp1 = lambda/(lambda + 1.0);
227: /*
228: Get a pointer to vector data.
229: - For default PETSc vectors, VecGetArray() returns a pointer to
230: the data array. Otherwise, the routine is implementation dependent.
231: - You MUST call VecRestoreArray() when you no longer need access to
232: the array.
233: */
234: DMDAVecGetArray(user->da,X,&x);
236: /*
237: Get local grid boundaries (for 3-dimensional DMDA):
238: xs, ys, zs - starting grid indices (no ghost points)
239: xm, ym, zm - widths of local grid (no ghost points)
241: */
242: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
244: /*
245: Compute initial guess over the locally owned part of the grid
246: */
247: for (k=zs; k<zs+zm; k++) {
248: tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
249: for (j=ys; j<ys+ym; j++) {
250: tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
251: for (i=xs; i<xs+xm; i++) {
252: if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
253: /* boundary conditions are all zero Dirichlet */
254: x[k][j][i] = 0.0;
255: } else {
256: x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
257: }
258: }
259: }
260: }
262: /*
263: Restore vector
264: */
265: DMDAVecRestoreArray(user->da,X,&x);
266: return(0);
267: }
268: /* ------------------------------------------------------------------- */
269: /*
270: FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch
272: Input Parameters:
273: . snes - the SNES context
274: . localX - input vector, this contains the ghosted region needed
275: . ptr - optional user-defined context, as set by SNESSetFunction()
277: Output Parameter:
278: . F - function vector, this does not contain a ghosted region
279: */
280: PetscErrorCode FormFunctionLocal(SNES snes,Vec localX,Vec F,void *ptr)
281: {
282: AppCtx *user = (AppCtx*)ptr;
284: PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
285: PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
286: PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
287: DM da;
290: SNESGetDM(snes,&da);
291: 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);
293: lambda = user->param;
294: hx = 1.0/(PetscReal)(Mx-1);
295: hy = 1.0/(PetscReal)(My-1);
296: hz = 1.0/(PetscReal)(Mz-1);
297: sc = hx*hy*hz*lambda;
298: hxhzdhy = hx*hz/hy;
299: hyhzdhx = hy*hz/hx;
300: hxhydhz = hx*hy/hz;
302: /*
303: Get pointers to vector data
304: */
305: DMDAVecGetArrayRead(da,localX,&x);
306: DMDAVecGetArray(da,F,&f);
308: /*
309: Get local grid boundaries
310: */
311: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
313: /*
314: Compute function over the locally owned part of the grid
315: */
316: for (k=zs; k<zs+zm; k++) {
317: for (j=ys; j<ys+ym; j++) {
318: for (i=xs; i<xs+xm; i++) {
319: if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
320: f[k][j][i] = x[k][j][i];
321: } else {
322: u = x[k][j][i];
323: u_east = x[k][j][i+1];
324: u_west = x[k][j][i-1];
325: u_north = x[k][j+1][i];
326: u_south = x[k][j-1][i];
327: u_up = x[k+1][j][i];
328: u_down = x[k-1][j][i];
329: u_xx = (-u_east + two*u - u_west)*hyhzdhx;
330: u_yy = (-u_north + two*u - u_south)*hxhzdhy;
331: u_zz = (-u_up + two*u - u_down)*hxhydhz;
332: f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
333: }
334: }
335: }
336: }
338: /*
339: Restore vectors
340: */
341: DMDAVecRestoreArrayRead(da,localX,&x);
342: DMDAVecRestoreArray(da,F,&f);
343: PetscLogFlops(11.0*ym*xm);
344: return(0);
345: }
346: /* ------------------------------------------------------------------- */
347: /*
348: FormFunction - Evaluates nonlinear function, F(x) on the entire domain
350: Input Parameters:
351: . snes - the SNES context
352: . X - input vector
353: . ptr - optional user-defined context, as set by SNESSetFunction()
355: Output Parameter:
356: . F - function vector
357: */
358: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
359: {
361: Vec localX;
362: DM da;
365: SNESGetDM(snes,&da);
366: DMGetLocalVector(da,&localX);
368: /*
369: Scatter ghost points to local vector,using the 2-step process
370: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
371: By placing code between these two statements, computations can be
372: done while messages are in transition.
373: */
374: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
375: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
377: FormFunctionLocal(snes,localX,F,ptr);
378: DMRestoreLocalVector(da,&localX);
379: return(0);
380: }
381: /* ------------------------------------------------------------------- */
382: /*
383: FormJacobian - Evaluates Jacobian matrix.
385: Input Parameters:
386: . snes - the SNES context
387: . x - input vector
388: . ptr - optional user-defined context, as set by SNESSetJacobian()
390: Output Parameters:
391: . A - Jacobian matrix
392: . B - optionally different preconditioning matrix
394: */
395: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
396: {
397: AppCtx *user = (AppCtx*)ptr; /* user-defined application context */
398: Vec localX;
400: PetscInt i,j,k,Mx,My,Mz;
401: MatStencil col[7],row;
402: PetscInt xs,ys,zs,xm,ym,zm;
403: PetscScalar lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;
404: DM da;
407: SNESGetDM(snes,&da);
408: DMGetLocalVector(da,&localX);
409: 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);
411: lambda = user->param;
412: hx = 1.0/(PetscReal)(Mx-1);
413: hy = 1.0/(PetscReal)(My-1);
414: hz = 1.0/(PetscReal)(Mz-1);
415: sc = hx*hy*hz*lambda;
416: hxhzdhy = hx*hz/hy;
417: hyhzdhx = hy*hz/hx;
418: hxhydhz = hx*hy/hz;
420: /*
421: Scatter ghost points to local vector, using the 2-step process
422: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
423: By placing code between these two statements, computations can be
424: done while messages are in transition.
425: */
426: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
427: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
429: /*
430: Get pointer to vector data
431: */
432: DMDAVecGetArrayRead(da,localX,&x);
434: /*
435: Get local grid boundaries
436: */
437: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
439: /*
440: Compute entries for the locally owned part of the Jacobian.
441: - Currently, all PETSc parallel matrix formats are partitioned by
442: contiguous chunks of rows across the processors.
443: - Each processor needs to insert only elements that it owns
444: locally (but any non-local elements will be sent to the
445: appropriate processor during matrix assembly).
446: - Here, we set all entries for a particular row at once.
447: - We can set matrix entries either using either
448: MatSetValuesLocal() or MatSetValues(), as discussed above.
449: */
450: for (k=zs; k<zs+zm; k++) {
451: for (j=ys; j<ys+ym; j++) {
452: for (i=xs; i<xs+xm; i++) {
453: row.k = k; row.j = j; row.i = i;
454: /* boundary points */
455: if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
456: v[0] = 1.0;
457: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
458: } else {
459: /* interior grid points */
460: v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j; col[0].i = i;
461: v[1] = -hxhzdhy; col[1].k=k; col[1].j=j-1;col[1].i = i;
462: v[2] = -hyhzdhx; col[2].k=k; col[2].j=j; col[2].i = i-1;
463: v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i;
464: v[4] = -hyhzdhx; col[4].k=k; col[4].j=j; col[4].i = i+1;
465: v[5] = -hxhzdhy; col[5].k=k; col[5].j=j+1;col[5].i = i;
466: v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j; col[6].i = i;
467: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
468: }
469: }
470: }
471: }
472: DMDAVecRestoreArrayRead(da,localX,&x);
473: DMRestoreLocalVector(da,&localX);
475: /*
476: Assemble matrix, using the 2-step process:
477: MatAssemblyBegin(), MatAssemblyEnd().
478: */
479: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
480: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
482: /*
483: Normally since the matrix has already been assembled above; this
484: would do nothing. But in the matrix free mode -snes_mf_operator
485: this tells the "matrix-free" matrix that a new linear system solve
486: is about to be done.
487: */
489: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
490: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
492: /*
493: Tell the matrix we will never add a new nonzero location to the
494: matrix. If we do, it will generate an error.
495: */
496: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
497: return(0);
498: }
500: /*TEST
502: test:
503: nsize: 4
504: args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
506: test:
507: suffix: 2
508: nsize: 4
509: args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
511: test:
512: suffix: 3
513: nsize: 4
514: args: -fdcoloring -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
516: test:
517: suffix: 3_ds
518: nsize: 4
519: args: -fdcoloring -fdcoloring_ds -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
521: test:
522: suffix: 4
523: nsize: 4
524: args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1
525: requires: !single
527: TEST*/