Actual source code: ex14.c

  1: /*$Id: ex14.c,v 1.33 2001/08/07 21:30:50 bsmith Exp $*/

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

  5: static char help[] = "Solves a nonlinear system in parallel with a user-defined Newton method.\n\
  6: Uses SLES to solve the linearized Newton sytems.  This solver\n\
  7: is a very simplistic inexact Newton method.  The intent of this code is to\n\
  8: demonstrate the repeated solution of linear sytems with the same nonzero pattern.\n\
  9: \n\
 10: This is NOT the recommended approach for solving nonlinear problems with PETSc!\n\
 11: We urge users to employ the SNES component for solving nonlinear problems whenever\n\
 12: possible, as it offers many advantages over coding nonlinear solvers independently.\n\
 13: \n\
 14: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
 15: domain, using distributed arrays (DAs) to partition the parallel grid.\n\
 16: The command line options include:\n\
 17:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
 18:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
 19:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
 20:   -my <yg>, where <yg> = number of grid points in the y-direction\n\
 21:   -Nx <npx>, where <npx> = number of processors in the x-direction\n\
 22:   -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";

 24: /*T
 25:    Concepts: SLES^writing a user-defined nonlinear solver (parallel Bratu example);
 26:    Concepts: DA^using distributed arrays;
 27:    Processors: n
 28: T*/

 30: /* ------------------------------------------------------------------------

 32:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 33:     the partial differential equation
 34:   
 35:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 36:   
 37:     with boundary conditions
 38:    
 39:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 40:   
 41:     A finite difference approximation with the usual 5-point stencil
 42:     is used to discretize the boundary value problem to obtain a nonlinear 
 43:     system of equations.

 45:     The SNES version of this problem is:  snes/examples/tutorials/ex5.c
 46:     We urge users to employ the SNES component for solving nonlinear
 47:     problems whenever possible, as it offers many advantages over coding 
 48:     nonlinear solvers independently.

 50:   ------------------------------------------------------------------------- */

 52: /* 
 53:    Include "petscda.h" so that we can use distributed arrays (DAs).
 54:    Include "petscsles.h" so that we can use SLES solvers.  Note that this
 55:    file automatically includes:
 56:      petsc.h       - base PETSc routines   petscvec.h - vectors
 57:      petscsys.h    - system routines       petscmat.h - matrices
 58:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 59:      petscviewer.h - viewers               petscpc.h  - preconditioners
 60: */
 61:  #include petscda.h
 62:  #include petscsles.h

 64: /* 
 65:    User-defined application context - contains data needed by the 
 66:    application-provided call-back routines, ComputeJacobian() and
 67:    ComputeFunction().
 68: */
 69: typedef struct {
 70:    PetscReal   param;          /* test problem parameter */
 71:    int         mx,my;          /* discretization in x,y directions */
 72:    Vec         localX,localF; /* ghosted local vector */
 73:    DA          da;             /* distributed array data structure */
 74:    int         rank;           /* processor rank */
 75: } AppCtx;

 77: /* 
 78:    User-defined routines
 79: */
 80: extern int ComputeFunction(AppCtx*,Vec,Vec),FormInitialGuess(AppCtx*,Vec);
 81: extern int ComputeJacobian(AppCtx*,Vec,Mat,MatStructure*);

 85: int main(int argc,char **argv)
 86: {
 87:   /* -------------- Data to define application problem ---------------- */
 88:   MPI_Comm  comm;                /* communicator */
 89:   SLES      sles;                /* linear solver */
 90:   Vec       X,Y,F;             /* solution, update, residual vectors */
 91:   Mat       J;                   /* Jacobian matrix */
 92:   AppCtx    user;                /* user-defined work context */
 93:   int       Nx,Ny;              /* number of preocessors in x- and y- directions */
 94:   int       size;                /* number of processors */
 95:   PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
 96:   int       m,N,ierr;

 98:   /* --------------- Data to define nonlinear solver -------------- */
 99:   PetscReal       rtol = 1.e-8;        /* relative convergence tolerance */
100:   PetscReal       xtol = 1.e-8;        /* step convergence tolerance */
101:   PetscReal       ttol;                /* convergence tolerance */
102:   PetscReal       fnorm,ynorm,xnorm; /* various vector norms */
103:   int          max_nonlin_its = 10; /* maximum number of iterations for nonlinear solver */
104:   int          max_functions = 50;  /* maximum number of function evaluations */
105:   int          lin_its;             /* number of linear solver iterations for each step */
106:   int          i;                   /* nonlinear solve iteration number */
107:   MatStructure mat_flag;        /* flag indicating structure of preconditioner matrix */
108:   PetscTruth   no_output;           /* flag indicating whether to surpress output */
109:   PetscScalar  mone = -1.0;

111:   PetscInitialize(&argc,&argv,(char *)0,help);
112:   comm = PETSC_COMM_WORLD;
113:   MPI_Comm_rank(comm,&user.rank);
114:   PetscOptionsHasName(PETSC_NULL,"-no_output",&no_output);

116:   /*
117:      Initialize problem parameters
118:   */
119:   user.mx = 4; user.my = 4; user.param = 6.0;
120:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
121:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
122:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
123:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
124:     SETERRQ(1,"Lambda is out of range");
125:   }
126:   N = user.mx*user.my;

128:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129:      Create linear solver context
130:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

132:   SLESCreate(comm,&sles);

134:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135:      Create vector data structures
136:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

138:   /*
139:      Create distributed array (DA) to manage parallel grid and vectors
140:   */
141:   MPI_Comm_size(comm,&size);
142:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
143:   PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
144:   PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
145:   if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE))
146:     SETERRQ(1,"Incompatible number of processors:  Nx * Ny != size");
147:   DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,user.mx,
148:                     user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da);

150:   /*
151:      Extract global and local vectors from DA; then duplicate for remaining
152:      vectors that are the same types
153:   */
154:   DACreateGlobalVector(user.da,&X);
155:   DACreateLocalVector(user.da,&user.localX);
156:   VecDuplicate(X,&F);
157:   VecDuplicate(X,&Y);
158:   VecDuplicate(user.localX,&user.localF);


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

172:      Note: Here we only approximately preallocate storage space for the
173:      Jacobian.  See the users manual for a discussion of better techniques
174:      for preallocating matrix memory.
175:   */
176:   if (size == 1) {
177:     MatCreateSeqAIJ(comm,N,N,5,PETSC_NULL,&J);
178:   } else {
179:     VecGetLocalSize(X,&m);
180:     MatCreateMPIAIJ(comm,m,m,N,N,5,PETSC_NULL,3,PETSC_NULL,&J);
181:   }

183:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184:      Customize linear solver; set runtime options
185:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

187:   /*
188:      Set runtime options (e.g.,-ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
189:   */
190:   SLESSetFromOptions(sles);

192:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193:      Evaluate initial guess
194:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

196:   FormInitialGuess(&user,X);
197:   ComputeFunction(&user,X,F);   /* Compute F(X)    */
198:   VecNorm(F,NORM_2,&fnorm);     /* fnorm = || F || */
199:   ttol = fnorm*rtol;
200:   if (!no_output) PetscPrintf(comm,"Initial function norm = %g\n",fnorm);

202:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203:      Solve nonlinear system with a user-defined method
204:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

206:   /* 
207:       This solver is a very simplistic inexact Newton method, with no
208:       no damping strategies or bells and whistles. The intent of this code
209:       is  merely to demonstrate the repeated solution with SLES of linear
210:       sytems with the same nonzero structure.

212:       This is NOT the recommended approach for solving nonlinear problems
213:       with PETSc!  We urge users to employ the SNES component for solving
214:       nonlinear problems whenever possible with application codes, as it
215:       offers many advantages over coding nonlinear solvers independently.
216:    */

218:   for (i=0; i<max_nonlin_its; i++) {

220:     /* 
221:         Compute the Jacobian matrix.  See the comments in this routine for
222:         important information about setting the flag mat_flag.
223:      */
224:     ComputeJacobian(&user,X,J,&mat_flag);

226:     /* 
227:         Solve J Y = F, where J is the Jacobian matrix.
228:           - First, set the SLES linear operators.  Here the matrix that
229:             defines the linear system also serves as the preconditioning
230:             matrix.
231:           - Then solve the Newton system.
232:      */
233:     SLESSetOperators(sles,J,J,mat_flag);
234:     SLESSolve(sles,F,Y,&lin_its);

236:     /* 
237:        Compute updated iterate
238:      */
239:     VecNorm(Y,NORM_2,&ynorm);       /* ynorm = || Y || */
240:     VecAYPX(&mone,X,Y);             /* Y <- X - Y      */
241:     VecCopy(Y,X);                   /* X <- Y          */
242:     VecNorm(X,NORM_2,&xnorm);       /* xnorm = || X || */
243:     if (!no_output) {
244:       PetscPrintf(comm,"   linear solve iterations = %d, xnorm=%g, ynorm=%g\n",lin_its,xnorm,ynorm);
245:     }

247:     /* 
248:        Evaluate new nonlinear function
249:      */
250:     ComputeFunction(&user,X,F);     /* Compute F(X)    */
251:     VecNorm(F,NORM_2,&fnorm);       /* fnorm = || F || */
252:     if (!no_output) {
253:       PetscPrintf(comm,"Iteration %d, function norm = %g\n",i+1,fnorm);
254:     }

256:     /*
257:        Test for convergence
258:      */
259:     if (fnorm <= ttol) {
260:       if (!no_output) {
261:          PetscPrintf(comm,"Converged due to function norm %g < %g (relative tolerance)\n",fnorm,ttol);
262:       }
263:       break;
264:     }
265:     if (ynorm < xtol*(xnorm)) {
266:       if (!no_output) {
267:          PetscPrintf(comm,"Converged due to small update length: %g < %g * %g\n",ynorm,xtol,xnorm);
268:       }
269:       break;
270:     }
271:     if (i > max_functions) {
272:       if (!no_output) {
273:         PetscPrintf(comm,"Exceeded maximum number of function evaluations: %d > %d\n",i,max_functions);
274:       }
275:       break;
276:     }
277:   }
278:   PetscPrintf(comm,"Number of Newton iterations = %d\n",i+1);

280:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281:      Free work space.  All PETSc objects should be destroyed when they
282:      are no longer needed.
283:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

285:   MatDestroy(J);           VecDestroy(Y);
286:   VecDestroy(user.localX); VecDestroy(X);
287:   VecDestroy(user.localF); VecDestroy(F);
288:   SLESDestroy(sles);  DADestroy(user.da);
289:   PetscFinalize();

291:   return 0;
292: }
293: /* ------------------------------------------------------------------- */
296: /* 
297:    FormInitialGuess - Forms initial approximation.

299:    Input Parameters:
300:    user - user-defined application context
301:    X - vector

303:    Output Parameter:
304:    X - vector
305:  */
306: int FormInitialGuess(AppCtx *user,Vec X)
307: {
308:   int          i,j,row,mx,my,ierr,xs,ys,xm,ym,gxm,gym,gxs,gys;
309:   PetscReal    one = 1.0,lambda,temp1,temp,hx,hy;
310:   PetscScalar  *x;
311:   Vec          localX = user->localX;

313:   mx = user->mx;            my = user->my;            lambda = user->param;
314:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
315:   temp1 = lambda/(lambda + one);

317:   /*
318:      Get a pointer to vector data.
319:        - For default PETSc vectors, VecGetArray() returns a pointer to
320:          the data array.  Otherwise, the routine is implementation dependent.
321:        - You MUST call VecRestoreArray() when you no longer need access to
322:          the array.
323:   */
324:   VecGetArray(localX,&x);

326:   /*
327:      Get local grid boundaries (for 2-dimensional DA):
328:        xs, ys   - starting grid indices (no ghost points)
329:        xm, ym   - widths of local grid (no ghost points)
330:        gxs, gys - starting grid indices (including ghost points)
331:        gxm, gym - widths of local grid (including ghost points)
332:   */
333:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
334:   DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

336:   /*
337:      Compute initial guess over the locally owned part of the grid
338:   */
339:   for (j=ys; j<ys+ym; j++) {
340:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
341:     for (i=xs; i<xs+xm; i++) {
342:       row = i - gxs + (j - gys)*gxm;
343:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
344:         x[row] = 0.0;
345:         continue;
346:       }
347:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
348:     }
349:   }

351:   /*
352:      Restore vector
353:   */
354:   VecRestoreArray(localX,&x);

356:   /*
357:      Insert values into global vector
358:   */
359:   DALocalToGlobal(user->da,localX,INSERT_VALUES,X);
360:   return 0;
361: }
362: /* ------------------------------------------------------------------- */
365: /* 
366:    ComputeFunction - Evaluates nonlinear function, F(x).

368:    Input Parameters:
369: .  X - input vector
370: .  user - user-defined application context

372:    Output Parameter:
373: .  F - function vector
374:  */
375: int ComputeFunction(AppCtx *user,Vec X,Vec F)
376: {
377:   int          ierr,i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
378:   PetscReal    two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
379:   PetscScalar  u,uxx,uyy,*x,*f;
380:   Vec          localX = user->localX,localF = user->localF;

382:   mx = user->mx;            my = user->my;            lambda = user->param;
383:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
384:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

386:   /*
387:      Scatter ghost points to local vector, using the 2-step process
388:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
389:      By placing code between these two statements, computations can be
390:      done while messages are in transition.
391:   */
392:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
393:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

395:   /*
396:      Get pointers to vector data
397:   */
398:   VecGetArray(localX,&x);
399:   VecGetArray(localF,&f);

401:   /*
402:      Get local grid boundaries
403:   */
404:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
405:   DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

407:   /*
408:      Compute function over the locally owned part of the grid
409:   */
410:   for (j=ys; j<ys+ym; j++) {
411:     row = (j - gys)*gxm + xs - gxs - 1;
412:     for (i=xs; i<xs+xm; i++) {
413:       row++;
414:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
415:         f[row] = x[row];
416:         continue;
417:       }
418:       u = x[row];
419:       uxx = (two*u - x[row-1] - x[row+1])*hydhx;
420:       uyy = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
421:       f[row] = uxx + uyy - sc*PetscExpScalar(u);
422:     }
423:   }

425:   /*
426:      Restore vectors
427:   */
428:   VecRestoreArray(localX,&x);
429:   VecRestoreArray(localF,&f);

431:   /*
432:      Insert values into global vector
433:   */
434:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
435:   PetscLogFlops(11*ym*xm);
436:   return 0;
437: }
438: /* ------------------------------------------------------------------- */
441: /*
442:    ComputeJacobian - Evaluates Jacobian matrix.

444:    Input Parameters:
445: .  x - input vector
446: .  user - user-defined application context

448:    Output Parameters:
449: .  jac - Jacobian matrix
450: .  flag - flag indicating matrix structure

452:    Notes:
453:    Due to grid point reordering with DAs, we must always work
454:    with the local grid points, and then transform them to the new
455:    global numbering with the "ltog" mapping (via DAGetGlobalIndices()).
456:    We cannot work directly with the global numbers for the original
457:    uniprocessor grid!
458: */
459: int ComputeJacobian(AppCtx *user,Vec X,Mat jac,MatStructure *flag)
460: {
461:   Vec     localX = user->localX;   /* local vector */
462:   int     *ltog;                   /* local-to-global mapping */
463:   int     ierr,i,j,row,mx,my,col[5];
464:   int     nloc,xs,ys,xm,ym,gxs,gys,gxm,gym,grow;
465:   PetscScalar  two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;

467:   mx = user->mx;            my = user->my;            lambda = user->param;
468:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
469:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

471:   /*
472:      Scatter ghost points to local vector, using the 2-step process
473:         DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
474:      By placing code between these two statements, computations can be
475:      done while messages are in transition.
476:   */
477:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
478:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

480:   /*
481:      Get pointer to vector data
482:   */
483:   VecGetArray(localX,&x);

485:   /*
486:      Get local grid boundaries
487:   */
488:   DAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
489:   DAGetGhostCorners(user->da,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL);

491:   /*
492:      Get the global node numbers for all local nodes, including ghost points
493:   */
494:   DAGetGlobalIndices(user->da,&nloc,&ltog);

496:   /* 
497:      Compute entries for the locally owned part of the Jacobian.
498:       - Currently, all PETSc parallel matrix formats are partitioned by
499:         contiguous chunks of rows across the processors. The "grow"
500:         parameter computed below specifies the global row number 
501:         corresponding to each local grid point.
502:       - Each processor needs to insert only elements that it owns
503:         locally (but any non-local elements will be sent to the
504:         appropriate processor during matrix assembly). 
505:       - Always specify global row and columns of matrix entries.
506:       - Here, we set all entries for a particular row at once.
507:   */
508:   for (j=ys; j<ys+ym; j++) {
509:     row = (j - gys)*gxm + xs - gxs - 1;
510:     for (i=xs; i<xs+xm; i++) {
511:       row++;
512:       grow = ltog[row];
513:       /* boundary points */
514:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
515:         MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
516:         continue;
517:       }
518:       /* interior grid points */
519:       v[0] = -hxdhy; col[0] = ltog[row - gxm];
520:       v[1] = -hydhx; col[1] = ltog[row - 1];
521:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
522:       v[3] = -hydhx; col[3] = ltog[row + 1];
523:       v[4] = -hxdhy; col[4] = ltog[row + gxm];
524:       MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
525:     }
526:   }

528:   /* 
529:      Assemble matrix, using the 2-step process:
530:        MatAssemblyBegin(), MatAssemblyEnd().
531:      By placing code between these two statements, computations can be
532:      done while messages are in transition.
533:   */
534:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
535:   VecRestoreArray(localX,&x);
536:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

538:   /*
539:      Set flag to indicate that the Jacobian matrix retains an identical
540:      nonzero structure throughout all nonlinear iterations (although the
541:      values of the entries change). Thus, we can save some work in setting
542:      up the preconditioner (e.g., no need to redo symbolic factorization for
543:      ILU/ICC preconditioners).
544:       - If the nonzero structure of the matrix is different during
545:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
546:         must be used instead.  If you are unsure whether the matrix
547:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
548:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
549:         believes your assertion and does not check the structure
550:         of the matrix.  If you erroneously claim that the structure
551:         is the same when it actually is not, the new preconditioner
552:         will not function correctly.  Thus, use this optimization
553:         feature with caution!
554:   */
555:   *flag = SAME_NONZERO_PATTERN;
556:   return 0;
557: }