Actual source code: ex14.c

  2: /* Program usage:  mpirun -np <procs> ex14 [-help] [all PETSc options] */

  4: static char help[] = "Bratu nonlinear PDE in 3d.\n\
  5: We solve the  Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\
  6: domain, using distributed arrays (DAs) to partition the parallel grid.\n\
  7: The command line options include:\n\
  8:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  9:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";

 11: /*T
 12:    Concepts: SNES^parallel Bratu example
 13:    Concepts: DA^using distributed arrays;
 14:    Processors: n
 15: T*/

 17: /* ------------------------------------------------------------------------

 19:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 20:     the partial differential equation
 21:   
 22:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 23:   
 24:     with boundary conditions
 25:    
 26:              u = 0  for  x = 0, x = 1, y = 0, y = 1, z = 0, z = 1
 27:   
 28:     A finite difference approximation with the usual 7-point stencil
 29:     is used to discretize the boundary value problem to obtain a nonlinear 
 30:     system of equations.


 33:   ------------------------------------------------------------------------- */

 35: /* 
 36:    Include "petscda.h" so that we can use distributed arrays (DAs).
 37:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 38:    file automatically includes:
 39:      petsc.h       - base PETSc routines   petscvec.h - vectors
 40:      petscsys.h    - system routines       petscmat.h - matrices
 41:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 42:      petscviewer.h - viewers               petscpc.h  - preconditioners
 43:      petscksp.h   - linear solvers
 44: */
 45:  #include petscda.h
 46:  #include petscsnes.h


 49: /* 
 50:    User-defined application context - contains data needed by the 
 51:    application-provided call-back routines, FormJacobian() and
 52:    FormFunction().
 53: */
 54: typedef struct {
 55:    PetscReal   param;          /* test problem parameter */
 56:    DA          da;             /* distributed array data structure */
 57: } AppCtx;

 59: /* 
 60:    User-defined routines
 61: */

 67: int main(int argc,char **argv)
 68: {
 69:   SNES                   snes;                 /* nonlinear solver */
 70:   Vec                    x,r;                  /* solution, residual vectors */
 71:   Mat                    J;                    /* Jacobian matrix */
 72:   AppCtx                 user;                 /* user-defined work context */
 73:   PetscInt               its;                  /* iterations for convergence */
 74:   PetscTruth             matrix_free,coloring;
 75:   PetscErrorCode         ierr;
 76:   PetscReal              bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm;
 77:   MatFDColoring          matfdcoloring;

 79:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 80:      Initialize program
 81:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 83:   PetscInitialize(&argc,&argv,(char *)0,help);

 85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 86:      Initialize problem parameters
 87:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 88:   user.param = 6.0;
 89:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
 90:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
 91:     SETERRQ(1,"Lambda is out of range");
 92:   }

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Create nonlinear solver context
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 97:   SNESCreate(PETSC_COMM_WORLD,&snes);

 99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:      Create distributed array (DA) to manage parallel grid and vectors
101:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102:   DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,4,4,4,PETSC_DECIDE,PETSC_DECIDE,
103:                     PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&user.da);

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

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

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:      Create matrix data structure; set Jacobian evaluation routine

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

129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   PetscOptionsHasName(PETSC_NULL,"-snes_mf",&matrix_free);
131:   PetscOptionsHasName(PETSC_NULL,"-fdcoloring",&coloring);
132:   if (!matrix_free) {
133:     if (coloring) {
134:       ISColoring    iscoloring;

136:       DAGetColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);
137:       DAGetMatrix(user.da,MATAIJ,&J);
138:       MatFDColoringCreate(J,iscoloring,&matfdcoloring);
139:       ISColoringDestroy(iscoloring);
140:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
141:       MatFDColoringSetFromOptions(matfdcoloring);
142:       SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
143:     } else {
144:       DAGetMatrix(user.da,MATAIJ,&J);
145:       SNESSetJacobian(snes,J,J,FormJacobian,&user);
146:     }
147:   }

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:      Customize nonlinear solver; set runtime options
151:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152:   SNESSetFromOptions(snes);

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

163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164:      Solve nonlinear system
165:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166:   SNESSolve(snes,PETSC_NULL,x);
167:   SNESGetIterationNumber(snes,&its);

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

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

181:   if (!matrix_free) {
182:     MatDestroy(J);
183:   }
184:   if (coloring) {
185:     MatFDColoringDestroy(matfdcoloring);
186:   }
187:   VecDestroy(x);
188:   VecDestroy(r);
189:   SNESDestroy(snes);
190:   DADestroy(user.da);
191:   PetscFinalize();

193:   return(0);
194: }
195: /* ------------------------------------------------------------------- */
198: /* 
199:    FormInitialGuess - Forms initial approximation.

201:    Input Parameters:
202:    user - user-defined application context
203:    X - vector

205:    Output Parameter:
206:    X - vector
207:  */
208: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
209: {
210:   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
212:   PetscReal      lambda,temp1,hx,hy,hz,tempk,tempj;
213:   PetscScalar    ***x;

216:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
217:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

219:   lambda = user->param;
220:   hx     = 1.0/(PetscReal)(Mx-1);
221:   hy     = 1.0/(PetscReal)(My-1);
222:   hz     = 1.0/(PetscReal)(Mz-1);
223:   temp1  = lambda/(lambda + 1.0);

225:   /*
226:      Get a pointer to vector data.
227:        - For default PETSc vectors, VecGetArray() returns a pointer to
228:          the data array.  Otherwise, the routine is implementation dependent.
229:        - You MUST call VecRestoreArray() when you no longer need access to
230:          the array.
231:   */
232:   DAVecGetArray(user->da,X,&x);

234:   /*
235:      Get local grid boundaries (for 3-dimensional DA):
236:        xs, ys, zs   - starting grid indices (no ghost points)
237:        xm, ym, zm   - widths of local grid (no ghost points)

239:   */
240:   DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

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

260:   /*
261:      Restore vector
262:   */
263:   DAVecRestoreArray(user->da,X,&x);
264:   return(0);
265: }
266: /* ------------------------------------------------------------------- */
269: /* 
270:    FormFunction - Evaluates nonlinear function, F(x).

272:    Input Parameters:
273: .  snes - the SNES context
274: .  X - input vector
275: .  ptr - optional user-defined context, as set by SNESSetFunction()

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

290:   DAGetLocalVector(user->da,&localX);
291:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
292:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

294:   lambda = user->param;
295:   hx     = 1.0/(PetscReal)(Mx-1);
296:   hy     = 1.0/(PetscReal)(My-1);
297:   hz     = 1.0/(PetscReal)(Mz-1);
298:   sc     = hx*hy*hz*lambda;
299:   hxhzdhy = hx*hz/hy;
300:   hyhzdhx = hy*hz/hx;
301:   hxhydhz = hx*hy/hz;

303:   /*
304:      Scatter ghost points to local vector,using the 2-step process
305:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
306:      By placing code between these two statements, computations can be
307:      done while messages are in transition.
308:   */
309:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
310:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

312:   /*
313:      Get pointers to vector data
314:   */
315:   DAVecGetArray(user->da,localX,&x);
316:   DAVecGetArray(user->da,F,&f);

318:   /*
319:      Get local grid boundaries
320:   */
321:   DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

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

348:   /*
349:      Restore vectors
350:   */
351:   DAVecRestoreArray(user->da,localX,&x);
352:   DAVecRestoreArray(user->da,F,&f);
353:   DARestoreLocalVector(user->da,&localX);
354:   PetscLogFlops(11*ym*xm);
355:   return(0);
356: }
357: /* ------------------------------------------------------------------- */
360: /*
361:    FormJacobian - Evaluates Jacobian matrix.

363:    Input Parameters:
364: .  snes - the SNES context
365: .  x - input vector
366: .  ptr - optional user-defined context, as set by SNESSetJacobian()

368:    Output Parameters:
369: .  A - Jacobian matrix
370: .  B - optionally different preconditioning matrix
371: .  flag - flag indicating matrix structure

373: */
374: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
375: {
376:   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
377:   Mat            jac = *B;                /* Jacobian matrix */
378:   Vec            localX;
380:   PetscInt       i,j,k,Mx,My,Mz;
381:   MatStencil     col[7],row;
382:   PetscInt       xs,ys,zs,xm,ym,zm;
383:   PetscScalar    lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;


387:   DAGetLocalVector(user->da,&localX);
388:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
389:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

391:   lambda = user->param;
392:   hx     = 1.0/(PetscReal)(Mx-1);
393:   hy     = 1.0/(PetscReal)(My-1);
394:   hz     = 1.0/(PetscReal)(Mz-1);
395:   sc     = hx*hy*hz*lambda;
396:   hxhzdhy = hx*hz/hy;
397:   hyhzdhx = hy*hz/hx;
398:   hxhydhz = hx*hy/hz;

400:   /*
401:      Scatter ghost points to local vector, using the 2-step process
402:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
403:      By placing code between these two statements, computations can be
404:      done while messages are in transition.
405:   */
406:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
407:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

409:   /*
410:      Get pointer to vector data
411:   */
412:   DAVecGetArray(user->da,localX,&x);

414:   /*
415:      Get local grid boundaries
416:   */
417:   DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

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

455:   /* 
456:      Assemble matrix, using the 2-step process:
457:        MatAssemblyBegin(), MatAssemblyEnd().
458:   */
459:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
460:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

462:   /*
463:      Normally since the matrix has already been assembled above; this
464:      would do nothing. But in the matrix free mode -snes_mf_operator
465:      this tells the "matrix-free" matrix that a new linear system solve
466:      is about to be done.
467:   */

469:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
470:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

472:   /*
473:      Set flag to indicate that the Jacobian matrix retains an identical
474:      nonzero structure throughout all nonlinear iterations (although the
475:      values of the entries change). Thus, we can save some work in setting
476:      up the preconditioner (e.g., no need to redo symbolic factorization for
477:      ILU/ICC preconditioners).
478:       - If the nonzero structure of the matrix is different during
479:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
480:         must be used instead.  If you are unsure whether the matrix
481:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
482:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
483:         believes your assertion and does not check the structure
484:         of the matrix.  If you erroneously claim that the structure
485:         is the same when it actually is not, the new preconditioner
486:         will not function correctly.  Thus, use this optimization
487:         feature with caution!
488:   */
489:   *flag = SAME_NONZERO_PATTERN;


492:   /*
493:      Tell the matrix we will never add a new nonzero location to the
494:      matrix. If we do, it will generate an error.
495:   */
496:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
497:   return(0);
498: }