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 suppress output */

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

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

110:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, NULL));
111:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, NULL));
112:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, NULL));
113:   PetscCheck(user.param < bratu_lambda_max && user.param > bratu_lambda_min, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Lambda is out of range");
114:   N = user.mx * user.my;

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Create linear solver context
118:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

120:   PetscCall(KSPCreate(comm, &ksp));

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:      Create vector data structures
124:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

126:   /*
127:      Create distributed array (DMDA) to manage parallel grid and vectors
128:   */
129:   PetscCallMPI(MPI_Comm_size(comm, &size));
130:   Nx = PETSC_DECIDE;
131:   Ny = PETSC_DECIDE;
132:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-Nx", &Nx, NULL));
133:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-Ny", &Ny, NULL));
134:   PetscCheck(Nx * Ny == size || (Nx == PETSC_DECIDE && Ny == PETSC_DECIDE), PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "Incompatible number of processors:  Nx * Ny != size");
135:   PetscCall(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));
136:   PetscCall(DMSetFromOptions(user.da));
137:   PetscCall(DMSetUp(user.da));

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

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

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

170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171:      Customize linear solver; set runtime options
172:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

174:   /*
175:      Set runtime options (e.g.,-ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
176:   */
177:   PetscCall(KSPSetFromOptions(ksp));

179:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180:      Evaluate initial guess
181:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

189:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190:      Solve nonlinear system with a user-defined method
191:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

205:   for (i = 0; i < max_nonlin_its; i++) {
206:     /*
207:         Compute the Jacobian matrix.
208:      */
209:     PetscCall(ComputeJacobian(&user, X, J));

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

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

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

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

256:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
257:      Free work space.  All PETSc objects should be destroyed when they
258:      are no longer needed.
259:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

261:   PetscCall(MatDestroy(&J));
262:   PetscCall(VecDestroy(&Y));
263:   PetscCall(VecDestroy(&user.localX));
264:   PetscCall(VecDestroy(&X));
265:   PetscCall(VecDestroy(&F));
266:   PetscCall(KSPDestroy(&ksp));
267:   PetscCall(DMDestroy(&user.da));
268:   PetscCall(PetscFinalize());
269:   return 0;
270: }
271: /* ------------------------------------------------------------------- */
272: /*
273:    FormInitialGuess - Forms initial approximation.

275:    Input Parameters:
276:    user - user-defined application context
277:    X - vector

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

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

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

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

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

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

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

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

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

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

371:   /*
372:      Get pointers to vector data
373:   */
374:   PetscCall(VecGetArray(localX, &x));
375:   PetscCall(VecGetArray(F, &f));

377:   /*
378:      Get local grid boundaries
379:   */
380:   PetscCall(DMDAGetCorners(user->da, &xs, &ys, NULL, &xm, &ym, NULL));
381:   PetscCall(DMDAGetGhostCorners(user->da, &gxs, &gys, NULL, &gxm, &gym, NULL));

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

401:   /*
402:      Restore vectors
403:   */
404:   PetscCall(VecRestoreArray(localX, &x));
405:   PetscCall(VecRestoreArray(F, &f));
406:   PetscCall(PetscLogFlops(11.0 * ym * xm));
407:   return PETSC_SUCCESS;
408: }
409: /* ------------------------------------------------------------------- */
410: /*
411:    ComputeJacobian - Evaluates Jacobian matrix.

413:    Input Parameters:
414: .  x - input vector
415: .  user - user-defined application context

417:    Output Parameters:
418: +  jac - Jacobian matrix
419: -  flag - flag indicating matrix structure

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

437:   mx     = user->mx;
438:   my     = user->my;
439:   lambda = user->param;
440:   hx     = one / (PetscReal)(mx - 1);
441:   hy     = one / (PetscReal)(my - 1);
442:   sc     = hx * hy;
443:   hxdhy  = hx / hy;
444:   hydhx  = hy / hx;

446:   /*
447:      Scatter ghost points to local vector, using the 2-step process
448:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
449:      By placing code between these two statements, computations can be
450:      done while messages are in transition.
451:   */
452:   PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX));
453:   PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX));

455:   /*
456:      Get pointer to vector data
457:   */
458:   PetscCall(VecGetArray(localX, &x));

460:   /*
461:      Get local grid boundaries
462:   */
463:   PetscCall(DMDAGetCorners(user->da, &xs, &ys, NULL, &xm, &ym, NULL));
464:   PetscCall(DMDAGetGhostCorners(user->da, &gxs, &gys, NULL, &gxm, &gym, NULL));

466:   /*
467:      Get the global node numbers for all local nodes, including ghost points
468:   */
469:   PetscCall(DMGetLocalToGlobalMapping(user->da, &ltogm));
470:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &ltog));

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

510:   /*
511:      Assemble matrix, using the 2-step process:
512:        MatAssemblyBegin(), MatAssemblyEnd().
513:      By placing code between these two statements, computations can be
514:      done while messages are in transition.
515:   */
516:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
517:   PetscCall(VecRestoreArray(localX, &x));
518:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

520:   return PETSC_SUCCESS;
521: }

523: /*TEST

525:     test:

527: TEST*/