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:  mpirun -np <procs> ex5 [-help] [all PETSc options]
32:      e.g.,
33:       ./ex5 -fd_jacobian -mat_fd_coloring_view_draw -draw_pause -1
34:       mpirun -np 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:      petsc.h       - base PETSc routines   petscvec.h - vectors
43:      petscsys.h    - system routines       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;
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);

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

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

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

130:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131:   /* J can be type of MATAIJ,MATBAIJ or MATSBAIJ */
132:   DAGetMatrix(user.da,MATAIJ,&J);
133:
134:   A    = J;
135:   PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian",&fd_jacobian,0);
136:   PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian_ghosted",&fd_jacobian_ghosted,0);
146:   }
147: #endif

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

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

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

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

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

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

206:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
207:      Solve nonlinear system
208:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
209:   SNESSolve(snes,PETSC_NULL,x);
210:   SNESGetIterationNumber(snes,&its);

212:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
214:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);

216:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217:      Free work space.  All PETSc objects should be destroyed when they
218:      are no longer needed.
219:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

221:   if (A != J) {
222:     MatDestroy(A);
223:   }
224:   MatDestroy(J);
225:   if (matfdcoloring) {
226:     MatFDColoringDestroy(matfdcoloring);
227:   }
228:   VecDestroy(x);
229:   VecDestroy(r);
230:   SNESDestroy(snes);
232:   PetscFinalize();

234:   return(0);
235: }
236: /* ------------------------------------------------------------------- */
239: /*
240:    FormInitialGuess - Forms initial approximation.

242:    Input Parameters:
243:    user - user-defined application context
244:    X - vector

246:    Output Parameter:
247:    X - vector
248:  */
249: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
250: {
251:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
253:   PetscReal      lambda,temp1,temp,hx,hy;
254:   PetscScalar    **x;

257:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
258:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

260:   lambda = user->param;
261:   hx     = 1.0/(PetscReal)(Mx-1);
262:   hy     = 1.0/(PetscReal)(My-1);
263:   temp1  = lambda/(lambda + 1.0);

265:   /*
266:      Get a pointer to vector data.
267:        - For default PETSc vectors, VecGetArray() returns a pointer to
268:          the data array.  Otherwise, the routine is implementation dependent.
269:        - You MUST call VecRestoreArray() when you no longer need access to
270:          the array.
271:   */
272:   DAVecGetArray(user->da,X,&x);

274:   /*
275:      Get local grid boundaries (for 2-dimensional DA):
276:        xs, ys   - starting grid indices (no ghost points)
277:        xm, ym   - widths of local grid (no ghost points)

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

282:   /*
283:      Compute initial guess over the locally owned part of the grid
284:   */
285:   for (j=ys; j<ys+ym; j++) {
286:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
287:     for (i=xs; i<xs+xm; i++) {
288:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
289:         /* boundary conditions are all zero Dirichlet */
290:         x[j][i] = 0.0;
291:       } else {
292:         x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
293:       }
294:     }
295:   }

297:   /*
298:      Restore vector
299:   */
300:   DAVecRestoreArray(user->da,X,&x);

302:   return(0);
303: }
304: /* ------------------------------------------------------------------- */
307: /*
308:    FormFunctionLocal - Evaluates nonlinear function, F(x).

312:  */
313: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
314: {
316:   PetscInt       i,j;
317:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
318:   PetscScalar    u,uxx,uyy;

322:   lambda = user->param;
323:   hx     = 1.0/(PetscReal)(info->mx-1);
324:   hy     = 1.0/(PetscReal)(info->my-1);
325:   sc     = hx*hy*lambda;
326:   hxdhy  = hx/hy;
327:   hydhx  = hy/hx;
328:   /*
329:      Compute function over the locally owned part of the grid
330:   */
331:   for (j=info->ys; j<info->ys+info->ym; j++) {
332:     for (i=info->xs; i<info->xs+info->xm; i++) {
333:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
334:         f[j][i] = x[j][i];
335:       } else {
336:         u       = x[j][i];
337:         uxx     = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
338:         uyy     = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
339:         f[j][i] = uxx + uyy - sc*PetscExpScalar(u);
340:       }
341:     }
342:   }

344:   PetscLogFlops(11*info->ym*info->xm);
345:   return(0);
346: }

350: /*
351:    FormJacobianLocal - Evaluates Jacobian matrix.
352: */
353: PetscErrorCode FormJacobianLocal(DALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
354: {
356:   PetscInt       i,j;
357:   MatStencil     col[5],row;
358:   PetscScalar    lambda,v[5],hx,hy,hxdhy,hydhx,sc;

361:   lambda = user->param;
362:   hx     = 1.0/(PetscReal)(info->mx-1);
363:   hy     = 1.0/(PetscReal)(info->my-1);
364:   sc     = hx*hy*lambda;
365:   hxdhy  = hx/hy;
366:   hydhx  = hy/hx;

369:   /*
370:      Compute entries for the locally owned part of the Jacobian.
371:       - Currently, all PETSc parallel matrix formats are partitioned by
372:         contiguous chunks of rows across the processors.
373:       - Each processor needs to insert only elements that it owns
374:         locally (but any non-local elements will be sent to the
375:         appropriate processor during matrix assembly).
376:       - Here, we set all entries for a particular row at once.
377:       - We can set matrix entries either using either
378:         MatSetValuesLocal() or MatSetValues(), as discussed above.
379:   */
380:   for (j=info->ys; j<info->ys+info->ym; j++) {
381:     for (i=info->xs; i<info->xs+info->xm; i++) {
382:       row.j = j; row.i = i;
383:       /* boundary points */
384:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
385:         v[0] = 1.0;
386:         MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
387:       } else {
388:       /* interior grid points */
389:         v[0] = -hxdhy;                                           col[0].j = j - 1; col[0].i = i;
390:         v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
391:         v[2] = 2.0*(hydhx + hxdhy) - sc*PetscExpScalar(x[j][i]); col[2].j = row.j; col[2].i = row.i;
392:         v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
393:         v[4] = -hxdhy;                                           col[4].j = j + 1; col[4].i = i;
394:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
395:       }
396:     }
397:   }

399:   /*
400:      Assemble matrix, using the 2-step process:
401:        MatAssemblyBegin(), MatAssemblyEnd().
402:   */
403:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
404:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
405:   /*
406:      Tell the matrix we will never add a new nonzero location to the
407:      matrix. If we do, it will generate an error.
408:   */
409:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
410:   return(0);
411: }

413: /*
414:       Variant of FormFunction() that computes the function in Matlab
415: */
416: #if defined(PETSC_HAVE_MATLAB_ENGINE)
417: PetscErrorCode FormFunctionMatlab(SNES snes,Vec X,Vec F,void *ptr)
418: {
419:   AppCtx         *user = (AppCtx*)ptr;
421:   PetscInt       Mx,My;
422:   PetscReal      lambda,hx,hy;
423:   Vec            localX,localF;
424:   MPI_Comm       comm;

427:   DAGetLocalVector(user->da,&localX);
428:   DAGetLocalVector(user->da,&localF);
429:   PetscObjectSetName((PetscObject)localX,"localX");
430:   PetscObjectSetName((PetscObject)localF,"localF");
431:   DAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
432:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

434:   lambda = user->param;
435:   hx     = 1.0/(PetscReal)(Mx-1);
436:   hy     = 1.0/(PetscReal)(My-1);

438:   PetscObjectGetComm((PetscObject)snes,&comm);
439:   /*
440:      Scatter ghost points to local vector,using the 2-step process
441:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
442:      By placing code between these two statements, computations can be
443:      done while messages are in transition.
444:   */
445:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
446:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
447:   PetscMatlabEnginePut(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localX);
448:   PetscMatlabEngineEvaluate(PETSC_MATLAB_ENGINE_(comm),"localF=ex5m(localX,%18.16e,%18.16e,%18.16e)",hx,hy,lambda);
449:   PetscMatlabEngineGet(PETSC_MATLAB_ENGINE_(comm),(PetscObject)localF);

451:   /*
452:      Insert values into global vector
453:   */
454:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
455:   DARestoreLocalVector(user->da,&localX);
456:   DARestoreLocalVector(user->da,&localF);
457:   return(0);
458: }
459: #endif

463: /*
464:   MySNESDefaultComputeJacobianColor - Computes the Jacobian using
465:     finite differences and coloring to exploit matrix sparsity.
466:     It is customized from SNESDefaultComputeJacobianColor.
467:     The input global vector x1 is scattered to x1_local
468:     which then is passed into MatFDColoringApply() for reducing the
469:     VecScatterBingin/End.
470: */
471: PetscErrorCode MySNESDefaultComputeJacobianColor(SNES snes,Vec x1,Mat *J,Mat *B,MatStructure *flag,void *ctx)
472: {
473:   MatFDColoring  color = (MatFDColoring) ctx;
475:   PetscInt       freq,it;
476:   Vec            f;
477:   PetscErrorCode (*ff)(void),(*fd)(void);
478:   void           *fctx;
479:   DA             da;
480:   Vec            x1_loc;

483:   MatFDColoringGetFrequency(color,&freq);
484:   SNESGetIterationNumber(snes,&it);

486:   if ((freq > 1) && ((it % freq))) {
487:     PetscInfo2(color,"Skipping Jacobian recomputation, it %D, freq %D\n",it,freq);
488:     *flag = SAME_PRECONDITIONER;
489:   } else {
490:     PetscInfo2(color,"Computing Jacobian, it %D, freq %D\n",it,freq);
491:     *flag = SAME_NONZERO_PATTERN;
492:     SNESGetFunction(snes,&f,(PetscErrorCode (**)(SNES,Vec,Vec,void*))&ff,0);
493:     MatFDColoringGetFunction(color,&fd,&fctx);
494:     if (fd == ff) { /* reuse function value computed in SNES */
495:       MatFDColoringSetF(color,f);
496:     }
497:     /* Now, get x1_loc and scatter global x1 onto x1_loc */
498:     da = *(DA*)fctx;
499:     DAGetLocalVector(da,&x1_loc);
500:     DAGlobalToLocalBegin(da,x1,INSERT_VALUES,x1_loc);
501:     DAGlobalToLocalEnd(da,x1,INSERT_VALUES,x1_loc);
502:     MatFDColoringApply(*B,color,x1_loc,flag,snes);
503:     DARestoreLocalVector(da,&x1_loc);
504:   }
505:   if (*J != *B) {
506:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
507:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
508:   }
509:   return(0);
510: }

```