Actual source code: ex5.c

  2: static char help[] = "Bratu nonlinear PDE in 2d.\n\
  3: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
  4: domain, using distributed arrays (DAs) 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: DA^using distributed arrays;
 12:    Concepts: IS coloirng types;
 13:    Processors: n
 14: T*/

 16: /* ------------------------------------------------------------------------

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

 31:     Program usage:  mpiexec -n <procs> ex5 [-help] [all PETSc options] 
 32:      e.g.,
 33:       ./ex5 -fd_jacobian -mat_fd_coloring_view_draw -draw_pause -1
 34:       mpiexec -n 2 ./ex5 -fd_jacobian_ghosted -log_summary

 36:   ------------------------------------------------------------------------- */

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

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

 61: /* 
 62:    User-defined routines
 63: */

 71: int main(int argc,char **argv)
 72: {
 73:   SNES                   snes;                 /* nonlinear solver */
 74:   Vec                    x,r;                  /* solution, residual vectors */
 75:   Mat                    A,J;                    /* Jacobian matrix */
 76:   AppCtx                 user;                 /* user-defined work context */
 77:   PetscInt               its;                  /* iterations for convergence */
 78:   PetscTruth             matlab_function = PETSC_FALSE;
 79:   PetscTruth             fd_jacobian = PETSC_FALSE,adic_jacobian=PETSC_FALSE,fd_jacobian_ghosted=PETSC_FALSE;
 80:   PetscTruth             adicmf_jacobian = PETSC_FALSE;
 81:   PetscErrorCode         ierr;
 82:   PetscReal              bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
 83:   MatFDColoring          matfdcoloring = 0;
 84:   ISColoring             iscoloring;

 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:      Initialize program
 88:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Initialize problem parameters
 94:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   user.param = 6.0;
 96:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
 97:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
 98:     SETERRQ(1,"Lambda is out of range");
 99:   }

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Create nonlinear solver context
103:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104:   SNESCreate(PETSC_COMM_WORLD,&snes);

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Create distributed array (DA) to manage parallel grid and vectors
108:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
110:                     1,1,PETSC_NULL,PETSC_NULL,&user.da);
111:   DASetUniformCoordinates(user.da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);

113:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Extract global vectors from DA; then duplicate for remaining
115:      vectors that are the same types
116:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   DACreateGlobalVector(user.da,&x);
118:   VecDuplicate(x,&r);

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Create matrix data structure; set Jacobian evaluation routine

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

131:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132:   /* J can be type of MATAIJ,MATBAIJ or MATSBAIJ */
133:   DAGetMatrix(user.da,MATAIJ,&J);
134: 
135:   A    = J;
136:   PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian",&fd_jacobian,0);
137:   PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian_ghosted",&fd_jacobian_ghosted,0);
138:   PetscOptionsGetTruth(PETSC_NULL,"-adic_jacobian",&adic_jacobian,0);
139:   PetscOptionsGetTruth(PETSC_NULL,"-adicmf_jacobian",&adicmf_jacobian,0);
140: #if defined(PETSC_HAVE_ADIC)
141:   if (adicmf_jacobian) {
142:     DASetLocalAdicMFFunction(user.da,admf_FormFunctionLocal);
143:     MatRegisterDAAD();
144:     MatCreateDAAD(user.da,&A);
145:     MatDAADSetSNES(A,snes);
146:     MatDAADSetCtx(A,&user);
147:   }
148: #endif

150:   if (fd_jacobian) {
151:     DAGetColoring(user.da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
152:     MatFDColoringCreate(J,iscoloring,&matfdcoloring);
153:     ISColoringDestroy(iscoloring);
154:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESDAFormFunction,&user);
155:     MatFDColoringSetFromOptions(matfdcoloring);
156:     SNESSetJacobian(snes,A,J,SNESDefaultComputeJacobianColor,matfdcoloring);
157:   } else if (fd_jacobian_ghosted) {
158:     DAGetColoring(user.da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
159:     MatFDColoringCreate(J,iscoloring,&matfdcoloring);
160:     ISColoringDestroy(iscoloring);
161:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESDAFormFunction,&user);
162:     MatFDColoringSetFromOptions(matfdcoloring);
163:     /* now, call a customized SNESDefaultComputeJacobianColor() */
164:     SNESSetJacobian(snes,A,J,MySNESDefaultComputeJacobianColor,matfdcoloring);
165: #if defined(PETSC_HAVE_ADIC)
166:   } else if (adic_jacobian) {
167:     DAGetColoring(user.da,IS_COLORING_GHOSTED,MATAIJ,&iscoloring);
168:     MatSetColoring(J,iscoloring);
169:     ISColoringDestroy(iscoloring);
170:     SNESSetJacobian(snes,A,J,SNESDAComputeJacobianWithAdic,&user);
171: #endif
172:   } else {
173:     SNESSetJacobian(snes,A,J,SNESDAComputeJacobian,&user);
174:   }

176:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177:      Set local function evaluation routine
178:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179:   DASetLocalFunction(user.da,(DALocalFunction1)FormFunctionLocal);
180:   DASetLocalJacobian(user.da,(DALocalFunction1)FormJacobianLocal);
181:   DASetLocalAdicFunction(user.da,ad_FormFunctionLocal);

183:   /* Decide which FormFunction to use */
184:   PetscOptionsGetTruth(PETSC_NULL,"-matlab_function",&matlab_function,0);

186:   SNESSetFunction(snes,r,SNESDAFormFunction,&user);
187: #if defined(PETSC_HAVE_MATLAB_ENGINE)
188:   if (matlab_function) {
189:     SNESSetFunction(snes,r,FormFunctionMatlab,&user);
190:   }
191: #endif

193:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194:      Customize nonlinear solver; set runtime options
195:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196:   SNESSetFromOptions(snes);

198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199:      Evaluate initial guess
200:      Note: The user should initialize the vector, x, with the initial guess
201:      for the nonlinear solver prior to calling SNESSolve().  In particular,
202:      to employ an initial guess of zero, the user should explicitly set
203:      this vector to zero by calling VecSet().
204:   */

206:   {
207:     PetscTruth test_appctx = PETSC_FALSE;
208:     PetscOptionsGetTruth(PETSC_NULL,"-test_appctx",&test_appctx,0);
209:     if (test_appctx) {
210:       AppCtx *puser;
211:       SNESSetApplicationContext(snes,&user);
212:       SNESGetApplicationContext(snes,(void **)&puser);
213:       FormInitialGuess(puser,x);
214:     } else {
215:       FormInitialGuess(&user,x);
216:     }
217:   }

219:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220:      Solve nonlinear system
221:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
222:   SNESSolve(snes,PETSC_NULL,x);
223:   SNESGetIterationNumber(snes,&its);

225:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
226:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
228:   {
229:     PetscViewer viewer;

231:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
232:     //PetscViewerSetType(viewer, PETSC_VIEWER_DRAW);
233:     //PetscViewerDrawSetInfo(viewer, PETSC_NULL, "Solution", PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
234:     PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
235:     PetscViewerFileSetName(viewer, "ex5_sol.vtk");
236:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
237:     DAView(user.da, viewer);
238:     VecView(x, viewer);
239:     PetscViewerDestroy(viewer);
240:   }

242:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243:      Free work space.  All PETSc objects should be destroyed when they
244:      are no longer needed.
245:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

247:   if (A != J) {
248:     MatDestroy(A);
249:   }
250:   MatDestroy(J);
251:   if (matfdcoloring) {
252:     MatFDColoringDestroy(matfdcoloring);
253:   }
254:   VecDestroy(x);
255:   VecDestroy(r);
256:   SNESDestroy(snes);
257:   DADestroy(user.da);
258:   PetscFinalize();

260:   return(0);
261: }
262: /* ------------------------------------------------------------------- */
265: /* 
266:    FormInitialGuess - Forms initial approximation.

268:    Input Parameters:
269:    user - user-defined application context
270:    X - vector

272:    Output Parameter:
273:    X - vector
274:  */
275: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
276: {
277:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
279:   PetscReal      lambda,temp1,temp,hx,hy;
280:   PetscScalar    **x;

283:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
284:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

286:   lambda = user->param;
287:   hx     = 1.0/(PetscReal)(Mx-1);
288:   hy     = 1.0/(PetscReal)(My-1);
289:   temp1  = lambda/(lambda + 1.0);

291:   /*
292:      Get a pointer to vector data.
293:        - For default PETSc vectors, VecGetArray() returns a pointer to
294:          the data array.  Otherwise, the routine is implementation dependent.
295:        - You MUST call VecRestoreArray() when you no longer need access to
296:          the array.
297:   */
298:   DAVecGetArray(user->da,X,&x);

300:   /*
301:      Get local grid boundaries (for 2-dimensional DA):
302:        xs, ys   - starting grid indices (no ghost points)
303:        xm, ym   - widths of local grid (no ghost points)

305:   */
306:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

308:   /*
309:      Compute initial guess over the locally owned part of the grid
310:   */
311:   for (j=ys; j<ys+ym; j++) {
312:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
313:     for (i=xs; i<xs+xm; i++) {
314:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
315:         /* boundary conditions are all zero Dirichlet */
316:         x[j][i] = 0.0;
317:       } else {
318:         x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
319:       }
320:     }
321:   }

323:   /*
324:      Restore vector
325:   */
326:   DAVecRestoreArray(user->da,X,&x);

328:   return(0);
329: }
330: /* ------------------------------------------------------------------- */
333: /* 
334:    FormFunctionLocal - Evaluates nonlinear function, F(x).

336:        Process adiC(36): FormFunctionLocal

338:  */
339: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
340: {
342:   PetscInt       i,j;
343:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
344:   PetscScalar    u,uxx,uyy;


348:   lambda = user->param;
349:   hx     = 1.0/(PetscReal)(info->mx-1);
350:   hy     = 1.0/(PetscReal)(info->my-1);
351:   sc     = hx*hy*lambda;
352:   hxdhy  = hx/hy;
353:   hydhx  = hy/hx;
354:   /*
355:      Compute function over the locally owned part of the grid
356:   */
357:   for (j=info->ys; j<info->ys+info->ym; j++) {
358:     for (i=info->xs; i<info->xs+info->xm; i++) {
359:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
360:         f[j][i] = x[j][i];
361:       } else {
362:         u       = x[j][i];
363:         uxx     = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
364:         uyy     = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
365:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
366:       }
367:     }
368:   }

370:   PetscLogFlops(11.0*info->ym*info->xm);
371:   return(0);
372: }

376: /*
377:    FormJacobianLocal - Evaluates Jacobian matrix.
378: */
379: PetscErrorCode FormJacobianLocal(DALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
380: {
382:   PetscInt       i,j;
383:   MatStencil     col[5],row;
384:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;

387:   lambda = user->param;
388:   hx     = 1.0/(PetscReal)(info->mx-1);
389:   hy     = 1.0/(PetscReal)(info->my-1);
390:   sc     = hx*hy*lambda;
391:   hxdhy  = hx/hy;
392:   hydhx  = hy/hx;


395:   /* 
396:      Compute entries for the locally owned part of the Jacobian.
397:       - Currently, all PETSc parallel matrix formats are partitioned by
398:         contiguous chunks of rows across the processors. 
399:       - Each processor needs to insert only elements that it owns
400:         locally (but any non-local elements will be sent to the
401:         appropriate processor during matrix assembly). 
402:       - Here, we set all entries for a particular row at once.
403:       - We can set matrix entries either using either
404:         MatSetValuesLocal() or MatSetValues(), as discussed above.
405:   */
406:   for (j=info->ys; j<info->ys+info->ym; j++) {
407:     for (i=info->xs; i<info->xs+info->xm; i++) {
408:       row.j = j; row.i = i;
409:       /* boundary points */
410:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
411:         v[0] = 1.0;
412:         MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
413:       } else {
414:       /* interior grid points */
415:         v[0] = -hxdhy;                                           col[0].j = j - 1; col[0].i = i;
416:         v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
417:         v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[2].j = row.j; col[2].i = row.i;
418:         v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
419:         v[4] = -hxdhy;                                           col[4].j = j + 1; col[4].i = i;
420:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
421:       }
422:     }
423:   }

425:   /* 
426:      Assemble matrix, using the 2-step process:
427:        MatAssemblyBegin(), MatAssemblyEnd().
428:   */
429:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
430:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
431:   /*
432:      Tell the matrix we will never add a new nonzero location to the
433:      matrix. If we do, it will generate an error.
434:   */
435:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
436:   return(0);
437: }

439: /*
440:       Variant of FormFunction() that computes the function in Matlab
441: */
442: #if defined(PETSC_HAVE_MATLAB_ENGINE)
443: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
444: {
445:   AppCtx         *user = (AppCtx*)ptr;
447:   PetscInt       Mx,My;
448:   PetscReal      lambda,hx,hy;
449:   Vec            localX,localF;
450:   MPI_Comm       comm;

453:   DAGetLocalVector(user->da,&localX);
454:   DAGetLocalVector(user->da,&localF);
455:   PetscObjectSetName((PetscObject)localX,"localX");
456:   PetscObjectSetName((PetscObject)localF,"localF");
457:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
458:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

460:   lambda = user->param;
461:   hx     = 1.0/(PetscReal)(Mx-1);
462:   hy     = 1.0/(PetscReal)(My-1);

464:   PetscObjectGetComm((PetscObject)snes,&comm);
465:   /*
466:      Scatter ghost points to local vector,using the 2-step process
467:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
468:      By placing code between these two statements, computations can be
469:      done while messages are in transition.
470:   */
471:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
472:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
473:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
474:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
475:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

477:   /*
478:      Insert values into global vector
479:   */
480:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
481:   DARestoreLocalVector(user->da,&localX);
482:   DARestoreLocalVector(user->da,&localF);
483:   return(0);
484: }
485: #endif

489: /*
490:   MySNESDefaultComputeJacobianColor - Computes the Jacobian using
491:     finite differences and coloring to exploit matrix sparsity. 
492:     It is customized from SNESDefaultComputeJacobianColor.
493:     The input global vector x1 is scattered to x1_local
494:     which then is passed into MatFDColoringApply() for reducing the
495:     VecScatterBingin/End.
496: */
497: PetscErrorCode MySNESDefaultComputeJacobianColor(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
498: {
499:   MatFDColoring  color = (MatFDColoring) ctx;
501:   Vec            f;
502:   PetscErrorCode (*ff)(void),(*fd)(void);
503:   void           *fctx;
504:   DA             da;
505:   Vec            x1_loc;

508:   *flag = SAME_NONZERO_PATTERN;
509:   SNESGetFunction(snes,&f,(PetscErrorCode (**)(SNES,Vec,Vec,void*))&ff,0);
510:   MatFDColoringGetFunction(color,&fd,&fctx);
511:   if (fd == ff) { /* reuse function value computed in SNES */
512:     MatFDColoringSetF(color,f);
513:   }
514:   /* Now, get x1_loc and scatter global x1 onto x1_loc */
515:   da = *(DA*)fctx;
516:   DAGetLocalVector(da,&x1_loc);
517:   DAGlobalToLocalBegin(da,x1,INSERT_VALUES,x1_loc);
518:   DAGlobalToLocalEnd(da,x1,INSERT_VALUES,x1_loc);
519:   MatFDColoringApply(*B,color,x1_loc,flag,snes);
520:   DARestoreLocalVector(da,&x1_loc);
521:   if (*J != *B) {
522:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
523:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
524:   }
525:   return(0);
526: }