Actual source code: ex14.c

petsc-master 2016-08-28
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*);

 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);

103:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104:      Extract global vectors from DMDA; then duplicate for remaining
105:      vectors that are the same types
106:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107:   DMCreateGlobalVector(user.da,&x);
108:   VecDuplicate(x,&r);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Set function evaluation routine and vector
112:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Create matrix data structure; set Jacobian evaluation routine

118:      Set Jacobian matrix data structure and default Jacobian evaluation
119:      routine. User can override with:
120:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
121:                 (unless user explicitly sets preconditioner)
122:      -snes_mf_operator : form preconditioning matrix as set by the user,
123:                          but use matrix-free approx for Jacobian-vector
124:                          products within Newton-Krylov method
125:      -fdcoloring : using finite differences with coloring to compute the Jacobian

127:      Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option
128:      below is to test the call to MatFDColoringSetType().
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
131:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL);
132:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL);
133:   PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL);
134:   if (!matrix_free) {
135:     DMSetMatType(user.da,MATAIJ);
136:     DMCreateMatrix(user.da,&J);
137:     if (coloring) {
138:       ISColoring iscoloring;
139:       if (!local_coloring) {
140:         DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);
141:         MatFDColoringCreate(J,iscoloring,&matfdcoloring);
142:         MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
143:       } else {
144:         DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring);
145:         MatFDColoringCreate(J,iscoloring,&matfdcoloring);
146:         MatFDColoringUseDM(J,matfdcoloring);
147:         MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunctionLocal,&user);
148:       }
149:       if (coloring_ds) {
150:         MatFDColoringSetType(matfdcoloring,MATMFFD_DS);
151:       }
152:       MatFDColoringSetFromOptions(matfdcoloring);
153:       MatFDColoringSetUp(J,iscoloring,matfdcoloring);
154:       SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
155:       ISColoringDestroy(&iscoloring);
156:     } else {
157:       SNESSetJacobian(snes,J,J,FormJacobian,&user);
158:     }
159:   }

161:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162:      Customize nonlinear solver; set runtime options
163:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164:   SNESSetDM(snes,user.da);
165:   SNESSetFromOptions(snes);

167:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168:      Evaluate initial guess
169:      Note: The user should initialize the vector, x, with the initial guess
170:      for the nonlinear solver prior to calling SNESSolve().  In particular,
171:      to employ an initial guess of zero, the user should explicitly set
172:      this vector to zero by calling VecSet().
173:   */
174:   FormInitialGuess(&user,x);

176:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177:      Solve nonlinear system
178:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179:   SNESSolve(snes,NULL,x);
180:   SNESGetIterationNumber(snes,&its);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Explicitly check norm of the residual of the solution
184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   FormFunction(snes,x,r,(void*)&user);
186:   VecNorm(r,NORM_2,&fnorm);
187:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm);

189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:      Free work space.  All PETSc objects should be destroyed when they
191:      are no longer needed.
192:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

194:   MatDestroy(&J);
195:   VecDestroy(&x);
196:   VecDestroy(&r);
197:   SNESDestroy(&snes);
198:   DMDestroy(&user.da);
199:   MatFDColoringDestroy(&matfdcoloring);
200:   PetscFinalize();
201:   return ierr;
202: }
203: /* ------------------------------------------------------------------- */
206: /*
207:    FormInitialGuess - Forms initial approximation.

209:    Input Parameters:
210:    user - user-defined application context
211:    X - vector

213:    Output Parameter:
214:    X - vector
215:  */
216: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
217: {
218:   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
220:   PetscReal      lambda,temp1,hx,hy,hz,tempk,tempj;
221:   PetscScalar    ***x;

224:   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);

226:   lambda = user->param;
227:   hx     = 1.0/(PetscReal)(Mx-1);
228:   hy     = 1.0/(PetscReal)(My-1);
229:   hz     = 1.0/(PetscReal)(Mz-1);
230:   temp1  = lambda/(lambda + 1.0);

232:   /*
233:      Get a pointer to vector data.
234:        - For default PETSc vectors, VecGetArray() returns a pointer to
235:          the data array.  Otherwise, the routine is implementation dependent.
236:        - You MUST call VecRestoreArray() when you no longer need access to
237:          the array.
238:   */
239:   DMDAVecGetArray(user->da,X,&x);

241:   /*
242:      Get local grid boundaries (for 3-dimensional DMDA):
243:        xs, ys, zs   - starting grid indices (no ghost points)
244:        xm, ym, zm   - widths of local grid (no ghost points)

246:   */
247:   DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

249:   /*
250:      Compute initial guess over the locally owned part of the grid
251:   */
252:   for (k=zs; k<zs+zm; k++) {
253:     tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
254:     for (j=ys; j<ys+ym; j++) {
255:       tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
256:       for (i=xs; i<xs+xm; i++) {
257:         if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
258:           /* boundary conditions are all zero Dirichlet */
259:           x[k][j][i] = 0.0;
260:         } else {
261:           x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
262:         }
263:       }
264:     }
265:   }

267:   /*
268:      Restore vector
269:   */
270:   DMDAVecRestoreArray(user->da,X,&x);
271:   return(0);
272: }
273: /* ------------------------------------------------------------------- */
276: /*
277:    FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch

279:    Input Parameters:
280: .  snes - the SNES context
281: .  localX - input vector, this contains the ghosted region needed 
282: .  ptr - optional user-defined context, as set by SNESSetFunction()

284:    Output Parameter:
285: .  F - function vector, this does not contain a ghosted region
286:  */
287: PetscErrorCode FormFunctionLocal(SNES snes,Vec localX,Vec F,void *ptr)
288: {
289:   AppCtx         *user = (AppCtx*)ptr;
291:   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
292:   PetscReal      two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
293:   PetscScalar    u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
294:   DM             da;

297:   SNESGetDM(snes,&da);
298:   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);

300:   lambda  = user->param;
301:   hx      = 1.0/(PetscReal)(Mx-1);
302:   hy      = 1.0/(PetscReal)(My-1);
303:   hz      = 1.0/(PetscReal)(Mz-1);
304:   sc      = hx*hy*hz*lambda;
305:   hxhzdhy = hx*hz/hy;
306:   hyhzdhx = hy*hz/hx;
307:   hxhydhz = hx*hy/hz;

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

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

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

345:   /*
346:      Restore vectors
347:   */
348:   DMDAVecRestoreArrayRead(da,localX,&x);
349:   DMDAVecRestoreArray(da,F,&f);
350:   PetscLogFlops(11.0*ym*xm);
351:   return(0);
352: }
353: /* ------------------------------------------------------------------- */
356: /*
357:    FormFunction - Evaluates nonlinear function, F(x) on the entire domain

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

364:    Output Parameter:
365: .  F - function vector
366:  */
367: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
368: {
370:   Vec            localX;
371:   DM             da;

374:   SNESGetDM(snes,&da);
375:   DMGetLocalVector(da,&localX);

377:   /*
378:      Scatter ghost points to local vector,using the 2-step process
379:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
380:      By placing code between these two statements, computations can be
381:      done while messages are in transition.
382:   */
383:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
384:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

386:   FormFunctionLocal(snes,localX,F,ptr);
387:   DMRestoreLocalVector(da,&localX);
388:   return(0);
389: }
390: /* ------------------------------------------------------------------- */
393: /*
394:    FormJacobian - Evaluates Jacobian matrix.

396:    Input Parameters:
397: .  snes - the SNES context
398: .  x - input vector
399: .  ptr - optional user-defined context, as set by SNESSetJacobian()

401:    Output Parameters:
402: .  A - Jacobian matrix
403: .  B - optionally different preconditioning matrix

405: */
406: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
407: {
408:   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
409:   Vec            localX;
411:   PetscInt       i,j,k,Mx,My,Mz;
412:   MatStencil     col[7],row;
413:   PetscInt       xs,ys,zs,xm,ym,zm;
414:   PetscScalar    lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;
415:   DM             da;

418:   SNESGetDM(snes,&da);
419:   DMGetLocalVector(da,&localX);
420:   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);

422:   lambda  = user->param;
423:   hx      = 1.0/(PetscReal)(Mx-1);
424:   hy      = 1.0/(PetscReal)(My-1);
425:   hz      = 1.0/(PetscReal)(Mz-1);
426:   sc      = hx*hy*hz*lambda;
427:   hxhzdhy = hx*hz/hy;
428:   hyhzdhx = hy*hz/hx;
429:   hxhydhz = hx*hy/hz;

431:   /*
432:      Scatter ghost points to local vector, using the 2-step process
433:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
434:      By placing code between these two statements, computations can be
435:      done while messages are in transition.
436:   */
437:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
438:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

440:   /*
441:      Get pointer to vector data
442:   */
443:   DMDAVecGetArrayRead(da,localX,&x);

445:   /*
446:      Get local grid boundaries
447:   */
448:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

450:   /*
451:      Compute entries for the locally owned part of the Jacobian.
452:       - Currently, all PETSc parallel matrix formats are partitioned by
453:         contiguous chunks of rows across the processors.
454:       - Each processor needs to insert only elements that it owns
455:         locally (but any non-local elements will be sent to the
456:         appropriate processor during matrix assembly).
457:       - Here, we set all entries for a particular row at once.
458:       - We can set matrix entries either using either
459:         MatSetValuesLocal() or MatSetValues(), as discussed above.
460:   */
461:   for (k=zs; k<zs+zm; k++) {
462:     for (j=ys; j<ys+ym; j++) {
463:       for (i=xs; i<xs+xm; i++) {
464:         row.k = k; row.j = j; row.i = i;
465:         /* boundary points */
466:         if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
467:           v[0] = 1.0;
468:           MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
469:         } else {
470:           /* interior grid points */
471:           v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j;  col[0].i = i;
472:           v[1] = -hxhzdhy; col[1].k=k;  col[1].j=j-1;col[1].i = i;
473:           v[2] = -hyhzdhx; col[2].k=k;  col[2].j=j;  col[2].i = i-1;
474:           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;
475:           v[4] = -hyhzdhx; col[4].k=k;  col[4].j=j;  col[4].i = i+1;
476:           v[5] = -hxhzdhy; col[5].k=k;  col[5].j=j+1;col[5].i = i;
477:           v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j;  col[6].i = i;
478:           MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
479:         }
480:       }
481:     }
482:   }
483:   DMDAVecRestoreArrayRead(da,localX,&x);
484:   DMRestoreLocalVector(da,&localX);

486:   /*
487:      Assemble matrix, using the 2-step process:
488:        MatAssemblyBegin(), MatAssemblyEnd().
489:   */
490:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
491:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

493:   /*
494:      Normally since the matrix has already been assembled above; this
495:      would do nothing. But in the matrix free mode -snes_mf_operator
496:      this tells the "matrix-free" matrix that a new linear system solve
497:      is about to be done.
498:   */

500:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
501:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

503:   /*
504:      Tell the matrix we will never add a new nonzero location to the
505:      matrix. If we do, it will generate an error.
506:   */
507:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
508:   return(0);
509: }