Actual source code: ex14.c

petsc-3.6.4 2016-04-12
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 FormFunction(SNES,Vec,Vec,void*),FormInitialGuess(AppCtx*,Vec);
 62: 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;                            /* Jacobian matrix */
 71:   AppCtx         user;                         /* user-defined work context */
 72:   PetscInt       its;                          /* iterations for convergence */
 73:   MatFDColoring  matfdcoloring;
 74:   PetscBool      matrix_free = PETSC_FALSE,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);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Initialize problem parameters
 86:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 87:   user.param = 6.0;
 88:   PetscOptionsGetReal(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,
100:                       PETSC_DECIDE,1,1,NULL,NULL,NULL,&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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127:   PetscOptionsGetBool(NULL,"-snes_mf",&matrix_free,NULL);
128:   PetscOptionsGetBool(NULL,"-fdcoloring",&coloring,NULL);
129:   if (!matrix_free) {
130:     DMSetMatType(user.da,MATAIJ);
131:     DMCreateMatrix(user.da,&J);
132:     if (coloring) {
133:       ISColoring iscoloring;
134:       DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);
135:       MatFDColoringCreate(J,iscoloring,&matfdcoloring);
136:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
137:       MatFDColoringSetFromOptions(matfdcoloring);
138:       MatFDColoringSetUp(J,iscoloring,matfdcoloring);
139:       SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
140:       ISColoringDestroy(&iscoloring);
141:     } else {
142:       SNESSetJacobian(snes,J,J,FormJacobian,&user);
143:     }
144:   }

146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147:      Customize nonlinear solver; set runtime options
148:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149:   SNESSetDM(snes,user.da);
150:   SNESSetFromOptions(snes);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:      Evaluate initial guess
154:      Note: The user should initialize the vector, x, with the initial guess
155:      for the nonlinear solver prior to calling SNESSolve().  In particular,
156:      to employ an initial guess of zero, the user should explicitly set
157:      this vector to zero by calling VecSet().
158:   */
159:   FormInitialGuess(&user,x);

161:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:      Solve nonlinear system
163:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164:   SNESSolve(snes,NULL,x);
165:   SNESGetIterationNumber(snes,&its);

167:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168:      Explicitly check norm of the residual of the solution
169:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170:   FormFunction(snes,x,r,(void*)&user);
171:   VecNorm(r,NORM_2,&fnorm);
172:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm);

174:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175:      Free work space.  All PETSc objects should be destroyed when they
176:      are no longer needed.
177:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

179:   if (!matrix_free) {
180:     MatDestroy(&J);
181:   }
182:   VecDestroy(&x);
183:   VecDestroy(&r);
184:   SNESDestroy(&snes);
185:   DMDestroy(&user.da);
186:   if (coloring) {MatFDColoringDestroy(&matfdcoloring);}
187:   PetscFinalize();
188:   return(0);
189: }
190: /* ------------------------------------------------------------------- */
193: /*
194:    FormInitialGuess - Forms initial approximation.

196:    Input Parameters:
197:    user - user-defined application context
198:    X - vector

200:    Output Parameter:
201:    X - vector
202:  */
203: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
204: {
205:   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: /* ------------------------------------------------------------------- */
263: /*
264:    FormFunction - Evaluates nonlinear function, F(x).

266:    Input Parameters:
267: .  snes - the SNES context
268: .  X - input vector
269: .  ptr - optional user-defined context, as set by SNESSetFunction()

271:    Output Parameter:
272: .  F - function vector
273:  */
274: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
275: {
276:   AppCtx         *user = (AppCtx*)ptr;
278:   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
279:   PetscReal      two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
280:   PetscScalar    u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
281:   Vec            localX;
282:   DM             da;

285:   SNESGetDM(snes,&da);
286:   DMGetLocalVector(da,&localX);
287:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
288:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

290:   lambda  = user->param;
291:   hx      = 1.0/(PetscReal)(Mx-1);
292:   hy      = 1.0/(PetscReal)(My-1);
293:   hz      = 1.0/(PetscReal)(Mz-1);
294:   sc      = hx*hy*hz*lambda;
295:   hxhzdhy = hx*hz/hy;
296:   hyhzdhx = hy*hz/hx;
297:   hxhydhz = hx*hy/hz;

299:   /*
300:      Scatter ghost points to local vector,using the 2-step process
301:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
302:      By placing code between these two statements, computations can be
303:      done while messages are in transition.
304:   */
305:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
306:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

308:   /*
309:      Get pointers to vector data
310:   */
311:   DMDAVecGetArrayRead(da,localX,&x);
312:   DMDAVecGetArray(da,F,&f);

314:   /*
315:      Get local grid boundaries
316:   */
317:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

319:   /*
320:      Compute function over the locally owned part of the grid
321:   */
322:   for (k=zs; k<zs+zm; k++) {
323:     for (j=ys; j<ys+ym; j++) {
324:       for (i=xs; i<xs+xm; i++) {
325:         if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
326:           f[k][j][i] = x[k][j][i];
327:         } else {
328:           u          = x[k][j][i];
329:           u_east     = x[k][j][i+1];
330:           u_west     = x[k][j][i-1];
331:           u_north    = x[k][j+1][i];
332:           u_south    = x[k][j-1][i];
333:           u_up       = x[k+1][j][i];
334:           u_down     = x[k-1][j][i];
335:           u_xx       = (-u_east + two*u - u_west)*hyhzdhx;
336:           u_yy       = (-u_north + two*u - u_south)*hxhzdhy;
337:           u_zz       = (-u_up + two*u - u_down)*hxhydhz;
338:           f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
339:         }
340:       }
341:     }
342:   }

344:   /*
345:      Restore vectors
346:   */
347:   DMDAVecRestoreArrayRead(da,localX,&x);
348:   DMDAVecRestoreArray(da,F,&f);
349:   DMRestoreLocalVector(da,&localX);
350:   PetscLogFlops(11.0*ym*xm);
351:   return(0);
352: }
353: /* ------------------------------------------------------------------- */
356: /*
357:    FormJacobian - Evaluates Jacobian matrix.

359:    Input Parameters:
360: .  snes - the SNES context
361: .  x - input vector
362: .  ptr - optional user-defined context, as set by SNESSetJacobian()

364:    Output Parameters:
365: .  A - Jacobian matrix
366: .  B - optionally different preconditioning matrix

368: */
369: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
370: {
371:   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
372:   Vec            localX;
374:   PetscInt       i,j,k,Mx,My,Mz;
375:   MatStencil     col[7],row;
376:   PetscInt       xs,ys,zs,xm,ym,zm;
377:   PetscScalar    lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;
378:   DM             da;

381:   SNESGetDM(snes,&da);
382:   DMGetLocalVector(da,&localX);
383:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
384:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

386:   lambda  = user->param;
387:   hx      = 1.0/(PetscReal)(Mx-1);
388:   hy      = 1.0/(PetscReal)(My-1);
389:   hz      = 1.0/(PetscReal)(Mz-1);
390:   sc      = hx*hy*hz*lambda;
391:   hxhzdhy = hx*hz/hy;
392:   hyhzdhx = hy*hz/hx;
393:   hxhydhz = hx*hy/hz;

395:   /*
396:      Scatter ghost points to local vector, using the 2-step process
397:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
398:      By placing code between these two statements, computations can be
399:      done while messages are in transition.
400:   */
401:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
402:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

404:   /*
405:      Get pointer to vector data
406:   */
407:   DMDAVecGetArrayRead(da,localX,&x);

409:   /*
410:      Get local grid boundaries
411:   */
412:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

414:   /*
415:      Compute entries for the locally owned part of the Jacobian.
416:       - Currently, all PETSc parallel matrix formats are partitioned by
417:         contiguous chunks of rows across the processors.
418:       - Each processor needs to insert only elements that it owns
419:         locally (but any non-local elements will be sent to the
420:         appropriate processor during matrix assembly).
421:       - Here, we set all entries for a particular row at once.
422:       - We can set matrix entries either using either
423:         MatSetValuesLocal() or MatSetValues(), as discussed above.
424:   */
425:   for (k=zs; k<zs+zm; k++) {
426:     for (j=ys; j<ys+ym; j++) {
427:       for (i=xs; i<xs+xm; i++) {
428:         row.k = k; row.j = j; row.i = i;
429:         /* boundary points */
430:         if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
431:           v[0] = 1.0;
432:           MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
433:         } else {
434:           /* interior grid points */
435:           v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j;  col[0].i = i;
436:           v[1] = -hxhzdhy; col[1].k=k;  col[1].j=j-1;col[1].i = i;
437:           v[2] = -hyhzdhx; col[2].k=k;  col[2].j=j;  col[2].i = i-1;
438:           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;
439:           v[4] = -hyhzdhx; col[4].k=k;  col[4].j=j;  col[4].i = i+1;
440:           v[5] = -hxhzdhy; col[5].k=k;  col[5].j=j+1;col[5].i = i;
441:           v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j;  col[6].i = i;
442:           MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
443:         }
444:       }
445:     }
446:   }
447:   DMDAVecRestoreArrayRead(da,localX,&x);
448:   DMRestoreLocalVector(da,&localX);

450:   /*
451:      Assemble matrix, using the 2-step process:
452:        MatAssemblyBegin(), MatAssemblyEnd().
453:   */
454:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
455:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

457:   /*
458:      Normally since the matrix has already been assembled above; this
459:      would do nothing. But in the matrix free mode -snes_mf_operator
460:      this tells the "matrix-free" matrix that a new linear system solve
461:      is about to be done.
462:   */

464:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
465:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

467:   /*
468:      Tell the matrix we will never add a new nonzero location to the
469:      matrix. If we do, it will generate an error.
470:   */
471:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
472:   return(0);
473: }