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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: PetscInitialize(&argc,&argv,(char*)0,help);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Initialize problem parameters
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: user.param = 6.0;
79: PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Create nonlinear solver context
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: SNESCreate(PETSC_COMM_WORLD,&snes);
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Create distributed array (DMDA) to manage parallel grid and vectors
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: 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);
91: DMSetFromOptions(user.da);
92: DMSetUp(user.da);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Extract global vectors from DMDA; then duplicate for remaining
95: vectors that are the same types
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: DMCreateGlobalVector(user.da,&x);
98: VecDuplicate(x,&r);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Set function evaluation routine and vector
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: SNESSetFunction(snes,r,FormFunction,(void*)&user);
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Create matrix data structure; set Jacobian evaluation routine
108: Set Jacobian matrix data structure and default Jacobian evaluation
109: routine. User can override with:
110: -snes_mf : matrix-free Newton-Krylov method with no preconditioning
111: (unless user explicitly sets preconditioner)
112: -snes_mf_operator : form preconditioning matrix as set by the user,
113: but use matrix-free approx for Jacobian-vector
114: products within Newton-Krylov method
115: -fdcoloring : using finite differences with coloring to compute the Jacobian
117: Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option
118: below is to test the call to MatFDColoringSetType().
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
121: PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL);
122: PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL);
123: PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL);
124: if (!matrix_free) {
125: DMSetMatType(user.da,MATAIJ);
126: DMCreateMatrix(user.da,&J);
127: if (coloring) {
128: ISColoring iscoloring;
129: if (!local_coloring) {
130: DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);
131: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
132: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
133: } else {
134: DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring);
135: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
136: MatFDColoringUseDM(J,matfdcoloring);
137: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunctionLocal,&user);
138: }
139: if (coloring_ds) {
140: MatFDColoringSetType(matfdcoloring,MATMFFD_DS);
141: }
142: MatFDColoringSetFromOptions(matfdcoloring);
143: MatFDColoringSetUp(J,iscoloring,matfdcoloring);
144: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
145: ISColoringDestroy(&iscoloring);
146: } else {
147: SNESSetJacobian(snes,J,J,FormJacobian,&user);
148: }
149: }
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Customize nonlinear solver; set runtime options
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: SNESSetDM(snes,user.da);
155: SNESSetFromOptions(snes);
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Evaluate initial guess
159: Note: The user should initialize the vector, x, with the initial guess
160: for the nonlinear solver prior to calling SNESSolve(). In particular,
161: to employ an initial guess of zero, the user should explicitly set
162: this vector to zero by calling VecSet().
163: */
164: FormInitialGuess(&user,x);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Solve nonlinear system
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: SNESSolve(snes,NULL,x);
170: SNESGetIterationNumber(snes,&its);
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Explicitly check norm of the residual of the solution
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: FormFunction(snes,x,r,(void*)&user);
176: VecNorm(r,NORM_2,&fnorm);
177: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm);
179: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180: Free work space. All PETSc objects should be destroyed when they
181: are no longer needed.
182: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184: MatDestroy(&J);
185: VecDestroy(&x);
186: VecDestroy(&r);
187: SNESDestroy(&snes);
188: DMDestroy(&user.da);
189: MatFDColoringDestroy(&matfdcoloring);
190: PetscFinalize();
191: return 0;
192: }
193: /* ------------------------------------------------------------------- */
194: /*
195: FormInitialGuess - Forms initial approximation.
197: Input Parameters:
198: user - user-defined application context
199: X - vector
201: Output Parameter:
202: X - vector
203: */
204: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
205: {
206: PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
207: PetscReal lambda,temp1,hx,hy,hz,tempk,tempj;
208: PetscScalar ***x;
211: 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);
213: lambda = user->param;
214: hx = 1.0/(PetscReal)(Mx-1);
215: hy = 1.0/(PetscReal)(My-1);
216: hz = 1.0/(PetscReal)(Mz-1);
217: temp1 = lambda/(lambda + 1.0);
219: /*
220: Get a pointer to vector data.
221: - For default PETSc vectors, VecGetArray() returns a pointer to
222: the data array. Otherwise, the routine is implementation dependent.
223: - You MUST call VecRestoreArray() when you no longer need access to
224: the array.
225: */
226: DMDAVecGetArray(user->da,X,&x);
228: /*
229: Get local grid boundaries (for 3-dimensional DMDA):
230: xs, ys, zs - starting grid indices (no ghost points)
231: xm, ym, zm - widths of local grid (no ghost points)
233: */
234: DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
236: /*
237: Compute initial guess over the locally owned part of the grid
238: */
239: for (k=zs; k<zs+zm; k++) {
240: tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
241: for (j=ys; j<ys+ym; j++) {
242: tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
243: for (i=xs; i<xs+xm; i++) {
244: if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
245: /* boundary conditions are all zero Dirichlet */
246: x[k][j][i] = 0.0;
247: } else {
248: x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
249: }
250: }
251: }
252: }
254: /*
255: Restore vector
256: */
257: DMDAVecRestoreArray(user->da,X,&x);
258: return 0;
259: }
260: /* ------------------------------------------------------------------- */
261: /*
262: FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch
264: Input Parameters:
265: . snes - the SNES context
266: . localX - input vector, this contains the ghosted region needed
267: . ptr - optional user-defined context, as set by SNESSetFunction()
269: Output Parameter:
270: . F - function vector, this does not contain a ghosted region
271: */
272: PetscErrorCode FormFunctionLocal(SNES snes,Vec localX,Vec F,void *ptr)
273: {
274: AppCtx *user = (AppCtx*)ptr;
275: PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
276: PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
277: PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
278: DM da;
281: SNESGetDM(snes,&da);
282: 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);
284: lambda = user->param;
285: hx = 1.0/(PetscReal)(Mx-1);
286: hy = 1.0/(PetscReal)(My-1);
287: hz = 1.0/(PetscReal)(Mz-1);
288: sc = hx*hy*hz*lambda;
289: hxhzdhy = hx*hz/hy;
290: hyhzdhx = hy*hz/hx;
291: hxhydhz = hx*hy/hz;
293: /*
294: Get pointers to vector data
295: */
296: DMDAVecGetArrayRead(da,localX,&x);
297: DMDAVecGetArray(da,F,&f);
299: /*
300: Get local grid boundaries
301: */
302: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
304: /*
305: Compute function over the locally owned part of the grid
306: */
307: for (k=zs; k<zs+zm; k++) {
308: for (j=ys; j<ys+ym; j++) {
309: for (i=xs; i<xs+xm; i++) {
310: if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
311: f[k][j][i] = x[k][j][i];
312: } else {
313: u = x[k][j][i];
314: u_east = x[k][j][i+1];
315: u_west = x[k][j][i-1];
316: u_north = x[k][j+1][i];
317: u_south = x[k][j-1][i];
318: u_up = x[k+1][j][i];
319: u_down = x[k-1][j][i];
320: u_xx = (-u_east + two*u - u_west)*hyhzdhx;
321: u_yy = (-u_north + two*u - u_south)*hxhzdhy;
322: u_zz = (-u_up + two*u - u_down)*hxhydhz;
323: f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
324: }
325: }
326: }
327: }
329: /*
330: Restore vectors
331: */
332: DMDAVecRestoreArrayRead(da,localX,&x);
333: DMDAVecRestoreArray(da,F,&f);
334: PetscLogFlops(11.0*ym*xm);
335: return 0;
336: }
337: /* ------------------------------------------------------------------- */
338: /*
339: FormFunction - Evaluates nonlinear function, F(x) on the entire domain
341: Input Parameters:
342: . snes - the SNES context
343: . X - input vector
344: . ptr - optional user-defined context, as set by SNESSetFunction()
346: Output Parameter:
347: . F - function vector
348: */
349: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
350: {
351: Vec localX;
352: DM da;
355: SNESGetDM(snes,&da);
356: DMGetLocalVector(da,&localX);
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(da,X,INSERT_VALUES,localX);
365: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
367: FormFunctionLocal(snes,localX,F,ptr);
368: DMRestoreLocalVector(da,&localX);
369: return 0;
370: }
371: /* ------------------------------------------------------------------- */
372: /*
373: FormJacobian - Evaluates Jacobian matrix.
375: Input Parameters:
376: . snes - the SNES context
377: . x - input vector
378: . ptr - optional user-defined context, as set by SNESSetJacobian()
380: Output Parameters:
381: . A - Jacobian matrix
382: . B - optionally different preconditioning matrix
384: */
385: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
386: {
387: AppCtx *user = (AppCtx*)ptr; /* user-defined application context */
388: Vec localX;
389: PetscInt i,j,k,Mx,My,Mz;
390: MatStencil col[7],row;
391: PetscInt xs,ys,zs,xm,ym,zm;
392: PetscScalar lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;
393: DM da;
396: SNESGetDM(snes,&da);
397: DMGetLocalVector(da,&localX);
398: 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);
400: lambda = user->param;
401: hx = 1.0/(PetscReal)(Mx-1);
402: hy = 1.0/(PetscReal)(My-1);
403: hz = 1.0/(PetscReal)(Mz-1);
404: sc = hx*hy*hz*lambda;
405: hxhzdhy = hx*hz/hy;
406: hyhzdhx = hy*hz/hx;
407: hxhydhz = hx*hy/hz;
409: /*
410: Scatter ghost points to local vector, using the 2-step process
411: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
412: By placing code between these two statements, computations can be
413: done while messages are in transition.
414: */
415: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
416: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
418: /*
419: Get pointer to vector data
420: */
421: DMDAVecGetArrayRead(da,localX,&x);
423: /*
424: Get local grid boundaries
425: */
426: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
428: /*
429: Compute entries for the locally owned part of the Jacobian.
430: - Currently, all PETSc parallel matrix formats are partitioned by
431: contiguous chunks of rows across the processors.
432: - Each processor needs to insert only elements that it owns
433: locally (but any non-local elements will be sent to the
434: appropriate processor during matrix assembly).
435: - Here, we set all entries for a particular row at once.
436: - We can set matrix entries either using either
437: MatSetValuesLocal() or MatSetValues(), as discussed above.
438: */
439: for (k=zs; k<zs+zm; k++) {
440: for (j=ys; j<ys+ym; j++) {
441: for (i=xs; i<xs+xm; i++) {
442: row.k = k; row.j = j; row.i = i;
443: /* boundary points */
444: if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
445: v[0] = 1.0;
446: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
447: } else {
448: /* interior grid points */
449: v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j; col[0].i = i;
450: v[1] = -hxhzdhy; col[1].k=k; col[1].j=j-1;col[1].i = i;
451: v[2] = -hyhzdhx; col[2].k=k; col[2].j=j; col[2].i = i-1;
452: 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;
453: v[4] = -hyhzdhx; col[4].k=k; col[4].j=j; col[4].i = i+1;
454: v[5] = -hxhzdhy; col[5].k=k; col[5].j=j+1;col[5].i = i;
455: v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j; col[6].i = i;
456: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
457: }
458: }
459: }
460: }
461: DMDAVecRestoreArrayRead(da,localX,&x);
462: DMRestoreLocalVector(da,&localX);
464: /*
465: Assemble matrix, using the 2-step process:
466: MatAssemblyBegin(), MatAssemblyEnd().
467: */
468: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
469: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
471: /*
472: Normally since the matrix has already been assembled above; this
473: would do nothing. But in the matrix free mode -snes_mf_operator
474: this tells the "matrix-free" matrix that a new linear system solve
475: is about to be done.
476: */
478: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
479: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
481: /*
482: Tell the matrix we will never add a new nonzero location to the
483: matrix. If we do, it will generate an error.
484: */
485: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
486: return 0;
487: }
489: /*TEST
491: test:
492: nsize: 4
493: args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
495: test:
496: suffix: 2
497: nsize: 4
498: args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
500: test:
501: suffix: 3
502: nsize: 4
503: args: -fdcoloring -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
505: test:
506: suffix: 3_ds
507: nsize: 4
508: args: -fdcoloring -fdcoloring_ds -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
510: test:
511: suffix: 4
512: nsize: 4
513: args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1
514: requires: !single
516: test:
517: suffix: 5
518: nsize: 4
519: args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1 -snes_type newtontrdc
520: requires: !single
522: TEST*/