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*/