Actual source code: ex14.c


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

 21: /* ------------------------------------------------------------------------

 23:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 24:     the partial differential equation

 26:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,

 28:     with boundary conditions

 30:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

 32:     A finite difference approximation with the usual 5-point stencil
 33:     is used to discretize the boundary value problem to obtain a nonlinear
 34:     system of equations.

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

 41:   ------------------------------------------------------------------------- */

 43: /*
 44:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 45:    Include "petscksp.h" so that we can use KSP solvers.  Note that this
 46:    file automatically includes:
 47:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 48:      petscmat.h - matrices
 49:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 50:      petscviewer.h - viewers               petscpc.h  - preconditioners
 51: */
 52: #include <petscdm.h>
 53: #include <petscdmda.h>
 54: #include <petscksp.h>

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

 68: /*
 69:    User-defined routines
 70: */
 71: extern PetscErrorCode ComputeFunction(AppCtx*,Vec,Vec),FormInitialGuess(AppCtx*,Vec);
 72: extern PetscErrorCode ComputeJacobian(AppCtx*,Vec,Mat);

 74: int main(int argc,char **argv)
 75: {
 76:   /* -------------- Data to define application problem ---------------- */
 77:   MPI_Comm       comm;                /* communicator */
 78:   KSP            ksp;                /* linear solver */
 79:   Vec            X,Y,F;             /* solution, update, residual vectors */
 80:   Mat            J;                   /* Jacobian matrix */
 81:   AppCtx         user;                /* user-defined work context */
 82:   PetscInt       Nx,Ny;              /* number of preocessors in x- and y- directions */
 83:   PetscMPIInt    size;                /* number of processors */
 84:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
 85:   PetscInt       m,N;

 87:   /* --------------- Data to define nonlinear solver -------------- */
 88:   PetscReal    rtol = 1.e-8;          /* relative convergence tolerance */
 89:   PetscReal    xtol = 1.e-8;          /* step convergence tolerance */
 90:   PetscReal    ttol;                  /* convergence tolerance */
 91:   PetscReal    fnorm,ynorm,xnorm;     /* various vector norms */
 92:   PetscInt     max_nonlin_its = 3;   /* maximum number of iterations for nonlinear solver */
 93:   PetscInt     max_functions  = 50;   /* maximum number of function evaluations */
 94:   PetscInt     lin_its;               /* number of linear solver iterations for each step */
 95:   PetscInt     i;                     /* nonlinear solve iteration number */
 96:   PetscBool    no_output = PETSC_FALSE;             /* flag indicating whether to surpress output */

 98:   PetscInitialize(&argc,&argv,(char*)0,help);
 99:   comm = PETSC_COMM_WORLD;
100:   PetscOptionsGetBool(NULL,NULL,"-no_output",&no_output,NULL);

102:   /*
103:      Initialize problem parameters
104:   */
105:   user.mx = 4; user.my = 4; user.param = 6.0;

107:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
108:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
109:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
111:   N = user.mx*user.my;

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114:      Create linear solver context
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

117:   KSPCreate(comm,&ksp);

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Create vector data structures
121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

123:   /*
124:      Create distributed array (DMDA) to manage parallel grid and vectors
125:   */
126:   MPI_Comm_size(comm,&size);
127:   Nx   = PETSC_DECIDE; Ny = PETSC_DECIDE;
128:   PetscOptionsGetInt(NULL,NULL,"-Nx",&Nx,NULL);
129:   PetscOptionsGetInt(NULL,NULL,"-Ny",&Ny,NULL);
131:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.da);
132:   DMSetFromOptions(user.da);
133:   DMSetUp(user.da);

135:   /*
136:      Extract global and local vectors from DMDA; then duplicate for remaining
137:      vectors that are the same types
138:   */
139:   DMCreateGlobalVector(user.da,&X);
140:   DMCreateLocalVector(user.da,&user.localX);
141:   VecDuplicate(X,&F);
142:   VecDuplicate(X,&Y);

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:      Create matrix data structure for Jacobian
146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147:   /*
148:      Note:  For the parallel case, vectors and matrices MUST be partitioned
149:      accordingly.  When using distributed arrays (DMDAs) to create vectors,
150:      the DMDAs determine the problem partitioning.  We must explicitly
151:      specify the local matrix dimensions upon its creation for compatibility
152:      with the vector distribution.  Thus, the generic MatCreate() routine
153:      is NOT sufficient when working with distributed arrays.

155:      Note: Here we only approximately preallocate storage space for the
156:      Jacobian.  See the users manual for a discussion of better techniques
157:      for preallocating matrix memory.
158:   */
159:   if (size == 1) {
160:     MatCreateSeqAIJ(comm,N,N,5,NULL,&J);
161:   } else {
162:     VecGetLocalSize(X,&m);
163:     MatCreateAIJ(comm,m,m,N,N,5,NULL,3,NULL,&J);
164:   }

166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:      Customize linear solver; set runtime options
168:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

170:   /*
171:      Set runtime options (e.g.,-ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
172:   */
173:   KSPSetFromOptions(ksp);

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:      Evaluate initial guess
177:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

179:   FormInitialGuess(&user,X);
180:   ComputeFunction(&user,X,F);   /* Compute F(X)    */
181:   VecNorm(F,NORM_2,&fnorm);     /* fnorm = || F || */
182:   ttol = fnorm*rtol;
183:   if (!no_output) PetscPrintf(comm,"Initial function norm = %g\n",(double)fnorm);

185:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186:      Solve nonlinear system with a user-defined method
187:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

189:   /*
190:       This solver is a very simplistic inexact Newton method, with no
191:       no damping strategies or bells and whistles. The intent of this code
192:       is  merely to demonstrate the repeated solution with KSP of linear
193:       systems with the same nonzero structure.

195:       This is NOT the recommended approach for solving nonlinear problems
196:       with PETSc!  We urge users to employ the SNES component for solving
197:       nonlinear problems whenever possible with application codes, as it
198:       offers many advantages over coding nonlinear solvers independently.
199:    */

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

203:     /*
204:         Compute the Jacobian matrix.
205:      */
206:     ComputeJacobian(&user,X,J);

208:     /*
209:         Solve J Y = F, where J is the Jacobian matrix.
210:           - First, set the KSP linear operators.  Here the matrix that
211:             defines the linear system also serves as the preconditioning
212:             matrix.
213:           - Then solve the Newton system.
214:      */
215:     KSPSetOperators(ksp,J,J);
216:     KSPSolve(ksp,F,Y);
217:     KSPGetIterationNumber(ksp,&lin_its);

219:     /*
220:        Compute updated iterate
221:      */
222:     VecNorm(Y,NORM_2,&ynorm);       /* ynorm = || Y || */
223:     VecAYPX(Y,-1.0,X);              /* Y <- X - Y      */
224:     VecCopy(Y,X);                   /* X <- Y          */
225:     VecNorm(X,NORM_2,&xnorm);       /* xnorm = || X || */
226:     if (!no_output) {
227:       PetscPrintf(comm,"   linear solve iterations = %D, xnorm=%g, ynorm=%g\n",lin_its,(double)xnorm,(double)ynorm);
228:     }

230:     /*
231:        Evaluate new nonlinear function
232:      */
233:     ComputeFunction(&user,X,F);     /* Compute F(X)    */
234:     VecNorm(F,NORM_2,&fnorm);       /* fnorm = || F || */
235:     if (!no_output) {
236:       PetscPrintf(comm,"Iteration %D, function norm = %g\n",i+1,(double)fnorm);
237:     }

239:     /*
240:        Test for convergence
241:      */
242:     if (fnorm <= ttol) {
243:       if (!no_output) {
244:         PetscPrintf(comm,"Converged due to function norm %g < %g (relative tolerance)\n",(double)fnorm,(double)ttol);
245:       }
246:       break;
247:     }
248:     if (ynorm < xtol*(xnorm)) {
249:       if (!no_output) {
250:         PetscPrintf(comm,"Converged due to small update length: %g < %g * %g\n",(double)ynorm,(double)xtol,(double)xnorm);
251:       }
252:       break;
253:     }
254:     if (i > max_functions) {
255:       if (!no_output) {
256:         PetscPrintf(comm,"Exceeded maximum number of function evaluations: %D > %D\n",i,max_functions);
257:       }
258:       break;
259:     }
260:   }
261:   PetscPrintf(comm,"Number of nonlinear iterations = %D\n",i);

263:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264:      Free work space.  All PETSc objects should be destroyed when they
265:      are no longer needed.
266:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

268:   MatDestroy(&J));           PetscCall(VecDestroy(&Y);
269:   VecDestroy(&user.localX)); PetscCall(VecDestroy(&X);
270:   VecDestroy(&F);
271:   KSPDestroy(&ksp));  PetscCall(DMDestroy(&user.da);
272:   PetscFinalize();
273:   return 0;
274: }
275: /* ------------------------------------------------------------------- */
276: /*
277:    FormInitialGuess - Forms initial approximation.

279:    Input Parameters:
280:    user - user-defined application context
281:    X - vector

283:    Output Parameter:
284:    X - vector
285:  */
286: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
287: {
288:   PetscInt    i,j,row,mx,my,xs,ys,xm,ym,gxm,gym,gxs,gys;
289:   PetscReal   one = 1.0,lambda,temp1,temp,hx,hy;
290:   PetscScalar *x;

292:   mx    = user->mx;            my = user->my;            lambda = user->param;
293:   hx    = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
294:   temp1 = lambda/(lambda + one);

296:   /*
297:      Get a pointer to vector data.
298:        - For default PETSc vectors, VecGetArray() returns a pointer to
299:          the data array.  Otherwise, the routine is implementation dependent.
300:        - You MUST call VecRestoreArray() when you no longer need access to
301:          the array.
302:   */
303:   VecGetArray(X,&x);

305:   /*
306:      Get local grid boundaries (for 2-dimensional DMDA):
307:        xs, ys   - starting grid indices (no ghost points)
308:        xm, ym   - widths of local grid (no ghost points)
309:        gxs, gys - starting grid indices (including ghost points)
310:        gxm, gym - widths of local grid (including ghost points)
311:   */
312:   DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
313:   DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);

315:   /*
316:      Compute initial guess over the locally owned part of the grid
317:   */
318:   for (j=ys; j<ys+ym; j++) {
319:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
320:     for (i=xs; i<xs+xm; i++) {
321:       row = i - gxs + (j - gys)*gxm;
322:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
323:         x[row] = 0.0;
324:         continue;
325:       }
326:       x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
327:     }
328:   }

330:   /*
331:      Restore vector
332:   */
333:   VecRestoreArray(X,&x);
334:   return 0;
335: }
336: /* ------------------------------------------------------------------- */
337: /*
338:    ComputeFunction - Evaluates nonlinear function, F(x).

340:    Input Parameters:
341: .  X - input vector
342: .  user - user-defined application context

344:    Output Parameter:
345: .  F - function vector
346:  */
347: PetscErrorCode ComputeFunction(AppCtx *user,Vec X,Vec F)
348: {
349:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,gxs,gys,gxm,gym;
350:   PetscReal      two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
351:   PetscScalar    u,uxx,uyy,*x,*f;
352:   Vec            localX = user->localX;

354:   mx = user->mx;            my = user->my;            lambda = user->param;
355:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
356:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

358:   /*
359:      Scatter ghost points to local vector, using the 2-step process
360:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
361:      By placing code between these two statements, computations can be
362:      done while messages are in transition.
363:   */
364:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
365:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

367:   /*
368:      Get pointers to vector data
369:   */
370:   VecGetArray(localX,&x);
371:   VecGetArray(F,&f);

373:   /*
374:      Get local grid boundaries
375:   */
376:   DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
377:   DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);

379:   /*
380:      Compute function over the locally owned part of the grid
381:   */
382:   for (j=ys; j<ys+ym; j++) {
383:     row = (j - gys)*gxm + xs - gxs - 1;
384:     for (i=xs; i<xs+xm; i++) {
385:       row++;
386:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
387:         f[row] = x[row];
388:         continue;
389:       }
390:       u      = x[row];
391:       uxx    = (two*u - x[row-1] - x[row+1])*hydhx;
392:       uyy    = (two*u - x[row-gxm] - x[row+gxm])*hxdhy;
393:       f[row] = uxx + uyy - sc*PetscExpScalar(u);
394:     }
395:   }

397:   /*
398:      Restore vectors
399:   */
400:   VecRestoreArray(localX,&x);
401:   VecRestoreArray(F,&f);
402:   PetscLogFlops(11.0*ym*xm);
403:   return 0;
404: }
405: /* ------------------------------------------------------------------- */
406: /*
407:    ComputeJacobian - Evaluates Jacobian matrix.

409:    Input Parameters:
410: .  x - input vector
411: .  user - user-defined application context

413:    Output Parameters:
414: .  jac - Jacobian matrix
415: .  flag - flag indicating matrix structure

417:    Notes:
418:    Due to grid point reordering with DMDAs, we must always work
419:    with the local grid points, and then transform them to the new
420:    global numbering with the "ltog" mapping
421:    We cannot work directly with the global numbers for the original
422:    uniprocessor grid!
423: */
424: PetscErrorCode ComputeJacobian(AppCtx *user,Vec X,Mat jac)
425: {
426:   Vec                    localX = user->localX;   /* local vector */
427:   const PetscInt         *ltog;                   /* local-to-global mapping */
428:   PetscInt               i,j,row,mx,my,col[5];
429:   PetscInt               xs,ys,xm,ym,gxs,gys,gxm,gym,grow;
430:   PetscScalar            two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;
431:   ISLocalToGlobalMapping ltogm;

433:   mx = user->mx;            my = user->my;            lambda = user->param;
434:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
435:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

437:   /*
438:      Scatter ghost points to local vector, using the 2-step process
439:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
440:      By placing code between these two statements, computations can be
441:      done while messages are in transition.
442:   */
443:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
444:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

446:   /*
447:      Get pointer to vector data
448:   */
449:   VecGetArray(localX,&x);

451:   /*
452:      Get local grid boundaries
453:   */
454:   DMDAGetCorners(user->da,&xs,&ys,NULL,&xm,&ym,NULL);
455:   DMDAGetGhostCorners(user->da,&gxs,&gys,NULL,&gxm,&gym,NULL);

457:   /*
458:      Get the global node numbers for all local nodes, including ghost points
459:   */
460:   DMGetLocalToGlobalMapping(user->da,&ltogm);
461:   ISLocalToGlobalMappingGetIndices(ltogm,&ltog);

463:   /*
464:      Compute entries for the locally owned part of the Jacobian.
465:       - Currently, all PETSc parallel matrix formats are partitioned by
466:         contiguous chunks of rows across the processors. The "grow"
467:         parameter computed below specifies the global row number
468:         corresponding to each local grid point.
469:       - Each processor needs to insert only elements that it owns
470:         locally (but any non-local elements will be sent to the
471:         appropriate processor during matrix assembly).
472:       - Always specify global row and columns of matrix entries.
473:       - Here, we set all entries for a particular row at once.
474:   */
475:   for (j=ys; j<ys+ym; j++) {
476:     row = (j - gys)*gxm + xs - gxs - 1;
477:     for (i=xs; i<xs+xm; i++) {
478:       row++;
479:       grow = ltog[row];
480:       /* boundary points */
481:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
482:         MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
483:         continue;
484:       }
485:       /* interior grid points */
486:       v[0] = -hxdhy; col[0] = ltog[row - gxm];
487:       v[1] = -hydhx; col[1] = ltog[row - 1];
488:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
489:       v[3] = -hydhx; col[3] = ltog[row + 1];
490:       v[4] = -hxdhy; col[4] = ltog[row + gxm];
491:       MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
492:     }
493:   }
494:   ISLocalToGlobalMappingRestoreIndices(ltogm,&ltog);

496:   /*
497:      Assemble matrix, using the 2-step process:
498:        MatAssemblyBegin(), MatAssemblyEnd().
499:      By placing code between these two statements, computations can be
500:      done while messages are in transition.
501:   */
502:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
503:   VecRestoreArray(localX,&x);
504:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

506:   return 0;
507: }

509: /*TEST

511:     test:

513: TEST*/