Actual source code: ex5.c

  2: static char help[] = "Solves a nonlinear system in parallel with SNES.\n\
  3: We solve the modified Bratu problem in a 2D rectangular domain,\n\
  4: using distributed arrays (DAs) to partition the parallel grid.\n\
  5: The command line options include:\n\
  6:   -lambda <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  7:   -kappa  <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  8:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
  9:   -my <yg>, where <yg> = number of grid points in the y-direction\n\
 10:   -Nx <npx>, where <npx> = number of processors in the x-direction\n\
 11:   -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";

 13: /*T
 14:    Concepts: SNES^solving a system of nonlinear equations (parallel Bratu example);
 15:    Concepts: DA^using distributed arrays;
 16:    Processors: n
 17: T*/

 19: /* ------------------------------------------------------------------------

 21:     Modified Solid Fuel Ignition problem.  This problem is modeled by
 22:     the partial differential equation

 24:         -Laplacian u - kappa*\PartDer{u}{x} - lambda*exp(u) = 0,

 26:     where

 28:          0 < x,y < 1,
 29:   
 30:     with boundary conditions
 31:    
 32:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 33:   
 34:     A finite difference approximation with the usual 5-point stencil
 35:     is used to discretize the boundary value problem to obtain a nonlinear 
 36:     system of equations.

 38:   ------------------------------------------------------------------------- */

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

 53: /* 
 54:    User-defined application context - contains data needed by the 
 55:    application-provided call-back routines, FormJacobian() and
 56:    FormFunction().
 57: */
 58: typedef struct {
 59:    PetscReal   param;          /* test problem parameter */
 60:    PetscReal   param2;         /* test problem parameter */
 61:    PetscInt    mx,my;          /* discretization in x, y directions */
 62:    Vec         localX,localF; /* ghosted local vector */
 63:    DA          da;             /* distributed array data structure */
 64:    PetscMPIInt rank;           /* processor rank */
 65: } AppCtx;

 67: /* 
 68:    User-defined routines
 69: */

 75: int main(int argc,char **argv)
 76: {
 77:   SNES           snes;                /* nonlinear solver */
 78:   Vec            x,r;                /* solution, residual vectors */
 79:   Mat            J;                   /* Jacobian matrix */
 80:   AppCtx         user;                /* user-defined work context */
 81:   PetscInt       its;                 /* iterations for convergence */
 82:   PetscInt       Nx,Ny;              /* number of preocessors in x- and y- directions */
 83:   PetscTruth     matrix_free = PETSC_FALSE;         /* flag - 1 indicates matrix-free version */
 84:   PetscMPIInt    size;                /* number of processors */
 85:   PetscInt       m,N;
 87:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
 88:   PetscReal      bratu_kappa_max = 10000,bratu_kappa_min = 0.;

 90:   PetscInitialize(&argc,&argv,(char *)0,help);
 91:   MPI_Comm_rank(PETSC_COMM_WORLD,&user.rank);

 93:   /*
 94:      Initialize problem parameters
 95:   */
 96:   user.mx = 4; user.my = 4; user.param = 6.0; user.param2 = 0.0;
 97:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
 98:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
 99:   PetscOptionsGetReal(PETSC_NULL,"-lambda",&user.param,PETSC_NULL);
100:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
101:     SETERRQ(1,"Lambda is out of range");
102:   }
103:   PetscOptionsGetReal(PETSC_NULL,"-kappa",&user.param2,PETSC_NULL);
104:   if (user.param2 >= bratu_kappa_max || user.param2 < bratu_kappa_min) {
105:     SETERRQ(1,"Kappa is out of range");
106:   }
107:   PetscPrintf(PETSC_COMM_WORLD,"Solving the Bratu problem with lambda=%G, kappa=%G\n",user.param,user.param2);

109:   N = user.mx*user.my;

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:      Create nonlinear solver context
113:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

115:   SNESCreate(PETSC_COMM_WORLD,&snes);

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:      Create vector data structures; set function evaluation routine
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   /*
122:      Create distributed array (DA) to manage parallel grid and vectors
123:   */
124:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
125:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
126:   PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
127:   PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
128:   if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE))
129:     SETERRQ(1,"Incompatible number of processors:  Nx * Ny != size");
130:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da);

132:   /*
133:      Visualize the distribution of the array across the processors
134:   */
135:   /*  DAView(user.da,PETSC_VIEWER_DRAW_WORLD); */


138:   /*
139:      Extract global and local vectors from DA; then duplicate for remaining
140:      vectors that are the same types
141:   */
142:   DACreateGlobalVector(user.da,&x);
143:   DACreateLocalVector(user.da,&user.localX);
144:   VecDuplicate(x,&r);
145:   VecDuplicate(user.localX,&user.localF);

147:   /* 
148:      Set function evaluation routine and vector
149:   */
150:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:      Create matrix data structure; set Jacobian evaluation routine
154:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

156:   /* 
157:      Set Jacobian matrix data structure and default Jacobian evaluation
158:      routine. User can override with:
159:      -snes_fd : default finite differencing approximation of Jacobian
160:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
161:                 (unless user explicitly sets preconditioner) 
162:      -snes_mf_operator : form preconditioning matrix as set by the user,
163:                          but use matrix-free approx for Jacobian-vector
164:                          products within Newton-Krylov method

166:      Note:  For the parallel case, vectors and matrices MUST be partitioned
167:      accordingly.  When using distributed arrays (DAs) to create vectors,
168:      the DAs determine the problem partitioning.  We must explicitly
169:      specify the local matrix dimensions upon its creation for compatibility
170:      with the vector distribution.  Thus, the generic MatCreate() routine
171:      is NOT sufficient when working with distributed arrays.

173:      Note: Here we only approximately preallocate storage space for the
174:      Jacobian.  See the users manual for a discussion of better techniques
175:      for preallocating matrix memory.
176:   */
177:   PetscOptionsGetTruth(PETSC_NULL,"-snes_mf",&matrix_free,PETSC_NULL);
178:   if (!matrix_free) {
179:     PetscTruth matrix_free_operator = PETSC_FALSE;
180:     PetscOptionsGetTruth(PETSC_NULL,"-snes_mf_operator",&matrix_free_operator,PETSC_NULL);
181:     if (matrix_free_operator) matrix_free = PETSC_FALSE;
182:   }
183:   if (!matrix_free) {
184:     if (size == 1) {
185:       MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,PETSC_NULL,&J);
186:     } else {
187:       VecGetLocalSize(x,&m);
188:       MatCreateMPIAIJ(PETSC_COMM_WORLD,m,m,N,N,5,PETSC_NULL,3,PETSC_NULL,&J);
189:     }
190:     SNESSetJacobian(snes,J,J,FormJacobian,&user);
191:   }

193:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194:      Customize nonlinear solver; set runtime options
195:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

197:   /*
198:      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
199:   */
200:   SNESSetFromOptions(snes);

202:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203:      Evaluate initial guess; then solve nonlinear system
204:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205:   /*
206:      Note: The user should initialize the vector, x, with the initial guess
207:      for the nonlinear solver prior to calling SNESSolve().  In particular,
208:      to employ an initial guess of zero, the user should explicitly set
209:      this vector to zero by calling VecSet().
210:   */
211:   FormInitialGuess(&user,x);
212:   SNESSolve(snes,PETSC_NULL,x);
213:   SNESGetIterationNumber(snes,&its);
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 (!matrix_free) {
222:     MatDestroy(J);
223:   }
224:   VecDestroy(user.localX); VecDestroy(x);
225:   VecDestroy(user.localF); VecDestroy(r);
226:   SNESDestroy(snes);  DADestroy(user.da);
227:   PetscFinalize();

229:   return 0;
230: }
231: /* ------------------------------------------------------------------- */
234: /* 
235:    FormInitialGuess - Forms initial approximation.

237:    Input Parameters:
238:    user - user-defined application context
239:    X - vector

241:    Output Parameter:
242:    X - vector
243:  */
244: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
245: {
246:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxm,gym,gxs,gys;
248:   PetscReal      one = 1.0,lambda,temp1,temp,hx,hy,hxdhy,hydhx,sc;
249:   PetscScalar    *x;
250:   Vec            localX = user->localX;

252:   mx = user->mx;            my = user->my;            lambda = user->param;
253:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
254:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;
255:   temp1 = lambda/(lambda + one);

257:   /*
258:      Get a pointer to vector data.
259:        - For default PETSc vectors,VecGetArray() returns a pointer to
260:          the data array.  Otherwise, the routine is implementation dependent.
261:        - You MUST call VecRestoreArray() when you no longer need access to
262:          the array.
263:   */
264:   VecGetArray(localX,&x);

266:   /*
267:      Get local grid boundaries (for 2-dimensional DA):
268:        xs, ys   - starting grid indices (no ghost points)
269:        xm, ym   - widths of local grid (no ghost points)
270:        gxs, gys - starting grid indices (including ghost points)
271:        gxm, gym - widths of local grid (including ghost points)
272:   */
273:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
274:   DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

276:   /*
277:      Compute initial guess over the locally owned part of the grid
278:   */
279:   for (j=ys; j<ys+ym; j++) {
280:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
281:     for (i=xs; i<xs+xm; i++) {
282:       row = i - gxs + (j - gys)*gxm;
283:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
284:         x[row] = 0.0;
285:         continue;
286:       }
287:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
288:     }
289:   }

291:   /*
292:      Restore vector
293:   */
294:   VecRestoreArray(localX,&x);

296:   /*
297:      Insert values into global vector
298:   */
299:   DALocalToGlobal(user->da,localX,INSERT_VALUES,X);
300:   return 0;
301: }
302: /* ------------------------------------------------------------------- */
305: /* 
306:    FormFunction - Evaluates nonlinear function, F(x).

308:    Input Parameters:
309: .  snes - the SNES context
310: .  X - input vector
311: .  ptr - optional user-defined context, as set by SNESSetFunction()

313:    Output Parameter:
314: .  F - function vector
315:  */
316: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
317: {
318:   AppCtx         *user = (AppCtx*)ptr;
320:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
321:   PetscReal      two = 2.0,one = 1.0,half = 0.5;
322:   PetscReal      lambda,hx,hy,hxdhy,hydhx,sc;
323:   PetscScalar    u,ux,uxx,uyy,*x,*f,kappa;
324:   Vec            localX = user->localX,localF = user->localF;

326:   mx = user->mx;            my = user->my;            lambda = user->param;
327:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
328:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;
329:   kappa = user->param2;

331:   /*
332:      Scatter ghost points to local vector, using the 2-step process
333:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
334:      By placing code between these two statements, computations can be
335:      done while messages are in transition.
336:   */
337:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
338:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

340:   /*
341:      Get pointers to vector data
342:   */
343:   VecGetArray(localX,&x);
344:   VecGetArray(localF,&f);

346:   /*
347:      Get local grid boundaries
348:   */
349:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
350:   DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

352:   /*
353:      Compute function over the locally owned part of the grid
354:   */
355:   for (j=ys; j<ys+ym; j++) {
356:     row = (j - gys)*gxm + xs - gxs - 1;
357:     for (i=xs; i<xs+xm; i++) {
358:       row++;
359:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
360:         f[row] = x[row];
361:         continue;
362:       }
363:       u = x[row];
364:       ux  = (x[row+1] - x[row-1])*half*hy;
365:       uxx = (two*u - x[row-1] - x[row+1])*hydhx;
366:       uyy = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
367:       f[row] = uxx + uyy - kappa*ux - sc*exp(u);
368:     }
369:   }

371:   /*
372:      Restore vectors
373:   */
374:   VecRestoreArray(localX,&x);
375:   VecRestoreArray(localF,&f);

377:   /*
378:      Insert values into global vector
379:   */
380:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
381:   PetscLogFlops(11.0*ym*xm);
382:   return 0;
383: }
384: /* ------------------------------------------------------------------- */
387: /*
388:    FormJacobian - Evaluates Jacobian matrix.

390:    Input Parameters:
391: .  snes - the SNES context
392: .  x - input vector
393: .  ptr - optional user-defined context, as set by SNESSetJacobian()

395:    Output Parameters:
396: .  A - Jacobian matrix
397: .  B - optionally different preconditioning matrix
398: .  flag - flag indicating matrix structure

400:    Notes:
401:    Due to grid point reordering with DAs, we must always work
402:    with the local grid points, and then transform them to the new
403:    global numbering with the "ltog" mapping (via DAGetGlobalIndices()).
404:    We cannot work directly with the global numbers for the original
405:    uniprocessor grid!
406: */
407: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
408: {
409:   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
410:   Mat            jac = *B;                /* Jacobian matrix */
411:   Vec            localX = user->localX;   /* local vector */
413:   PetscInt       *ltog;                   /* local-to-global mapping */
414:   PetscInt       i,j,row,mx,my,col[5];
415:   PetscInt       nloc,xs,ys,xm,ym,gxs,gys,gxm,gym,grow;
416:   PetscScalar    two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;

418:   mx = user->mx;            my = user->my;            lambda = user->param;
419:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
420:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

422:   /*
423:      Scatter ghost points to local vector,using the 2-step process
424:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
425:      By placing code between these two statements, computations can be
426:      done while messages are in transition.
427:   */
428:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
429:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

431:   /*
432:      Get pointer to vector data
433:   */
434:   VecGetArray(localX,&x);

436:   /*
437:      Get local grid boundaries
438:   */
439:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
440:   DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

442:   /*
443:      Get the global node numbers for all local nodes, including ghost points
444:   */
445:   DAGetGlobalIndices(user->da,&nloc,&ltog);

447:   /* 
448:      Compute entries for the locally owned part of the Jacobian.
449:       - Currently, all PETSc parallel matrix formats are partitioned by
450:         contiguous chunks of rows across the processors. The "grow"
451:         parameter computed below specifies the global row number 
452:         corresponding to each local grid point.
453:       - Each processor needs to insert only elements that it owns
454:         locally (but any non-local elements will be sent to the
455:         appropriate processor during matrix assembly). 
456:       - Always specify global row and columns of matrix entries.
457:       - Here, we set all entries for a particular row at once.
458:   */
459:   for (j=ys; j<ys+ym; j++) {
460:     row = (j - gys)*gxm + xs - gxs - 1;
461:     for (i=xs; i<xs+xm; i++) {
462:       row++;
463:       grow = ltog[row];
464:       /* boundary points */
465:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
466:         MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
467:         continue;
468:       }
469:       /* interior grid points */
470:       v[0] = -hxdhy; col[0] = ltog[row - gxm];
471:       v[1] = -hydhx; col[1] = ltog[row - 1];
472:       v[2] = two*(hydhx + hxdhy) - sc*lambda*exp(x[row]); col[2] = grow;
473:       v[3] = -hydhx; col[3] = ltog[row + 1];
474:       v[4] = -hxdhy; col[4] = ltog[row + gxm];
475:       MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
476:     }
477:   }

479:   /* 
480:      Assemble matrix, using the 2-step process:
481:        MatAssemblyBegin(), MatAssemblyEnd().
482:      By placing code between these two statements, computations can be
483:      done while messages are in transition.
484:   */
485:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
486:   VecRestoreArray(localX,&x);
487:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

489:   /*
490:      Set flag to indicate that the Jacobian matrix retains an identical
491:      nonzero structure throughout all nonlinear iterations (although the
492:      values of the entries change). Thus, we can save some work in setting
493:      up the preconditioner (e.g., no need to redo symbolic factorization for
494:      ILU/ICC preconditioners).
495:       - If the nonzero structure of the matrix is different during
496:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
497:         must be used instead.  If you are unsure whether the matrix
498:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
499:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
500:         believes your assertion and does not check the structure
501:         of the matrix.  If you erroneously claim that the structure
502:         is the same when it actually is not, the new preconditioner
503:         will not function correctly.  Thus, use this optimization
504:         feature with caution!
505:   */
506:   *flag = SAME_NONZERO_PATTERN;
507:   /*
508:       Tell the matrix we will never add a new nonzero location to the
509:     matrix. If we do it will generate an error.
510:   */
511:   /* MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); */
512:   return 0;
513: }