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