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