Actual source code: ex14.c

petsc-3.8.3 2017-12-09
Report Typos and Errors
```
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:   */
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:   */
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:   */

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:   }
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: }

```