Actual source code: eptorsion1.c

  1: /* Program usage: mpiexec -n 1 eptorsion1 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------

  5:   Elastic-plastic torsion problem.

  7:   The elastic plastic torsion problem arises from the determination
  8:   of the stress field on an infinitely long cylindrical bar, which is
  9:   equivalent to the solution of the following problem:

 11:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}

 13:   where C is the torsion angle per unit length.

 15:   The multiprocessor version of this code is eptorsion2.c.

 17: ---------------------------------------------------------------------- */

 19: /*
 20:   Include "petsctao.h" so that we can use TAO solvers.  Note that this
 21:   file automatically includes files for lower-level support, such as those
 22:   provided by the PETSc library:
 23:      petsc.h       - base PETSc routines   petscvec.h - vectors
 24:      petscsys.h    - system routines        petscmat.h - matrices
 25:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 26:      petscviewer.h - viewers               petscpc.h  - preconditioners
 27: */

 29: #include <petsctao.h>

 31: static char help[] = "Demonstrates use of the TAO package to solve \n\
 32: unconstrained minimization problems on a single processor.  This example \n\
 33: is based on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 \n\
 34: test suite.\n\
 35: The command line options are:\n\
 36:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 37:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 38:   -par <param>, where <param> = angle of twist per unit length\n\n";

 40: /*
 41:    User-defined application context - contains data needed by the
 42:    application-provided call-back routines, FormFunction(),
 43:    FormGradient(), and FormHessian().
 44: */

 46: typedef struct {
 47:   PetscReal param;      /* nonlinearity parameter */
 48:   PetscInt  mx, my;     /* discretization in x- and y-directions */
 49:   PetscInt  ndim;       /* problem dimension */
 50:   Vec       s, y, xvec; /* work space for computing Hessian */
 51:   PetscReal hx, hy;     /* mesh spacing in x- and y-directions */
 52: } AppCtx;

 54: /* -------- User-defined Routines --------- */

 56: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
 57: PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
 58: PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
 59: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
 60: PetscErrorCode HessianProductMat(Mat, Vec, Vec);
 61: PetscErrorCode HessianProduct(void *, Vec, Vec);
 62: PetscErrorCode MatrixFreeHessian(Tao, Vec, Mat, Mat, void *);
 63: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);

 65: PetscErrorCode main(int argc, char **argv)
 66: {
 67:   PetscInt    mx = 10; /* discretization in x-direction */
 68:   PetscInt    my = 10; /* discretization in y-direction */
 69:   Vec         x;       /* solution, gradient vectors */
 70:   PetscBool   flg;     /* A return value when checking for use options */
 71:   Tao         tao;     /* Tao solver context */
 72:   Mat         H;       /* Hessian matrix */
 73:   AppCtx      user;    /* application context */
 74:   PetscMPIInt size;    /* number of processes */
 75:   PetscReal   one = 1.0;

 77:   PetscBool test_lmvm = PETSC_FALSE;
 78:   KSP       ksp;
 79:   PC        pc;
 80:   Mat       M;
 81:   Vec       in, out, out2;
 82:   PetscReal mult_solve_dist;

 84:   /* Initialize TAO,PETSc */
 86:   PetscInitialize(&argc, &argv, (char *)0, help);
 87:   MPI_Comm_size(MPI_COMM_WORLD, &size);

 90:   /* Specify default parameters for the problem, check for command-line overrides */
 91:   user.param = 5.0;
 92:   PetscOptionsGetInt(NULL, NULL, "-my", &my, &flg);
 93:   PetscOptionsGetInt(NULL, NULL, "-mx", &mx, &flg);
 94:   PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg);
 95:   PetscOptionsGetBool(NULL, NULL, "-test_lmvm", &test_lmvm, &flg);

 97:   PetscPrintf(PETSC_COMM_SELF, "\n---- Elastic-Plastic Torsion Problem -----\n");
 98:   PetscPrintf(PETSC_COMM_SELF, "mx: %" PetscInt_FMT "     my: %" PetscInt_FMT "   \n\n", mx, my);
 99:   user.ndim = mx * my;
100:   user.mx   = mx;
101:   user.my   = my;
102:   user.hx   = one / (mx + 1);
103:   user.hy   = one / (my + 1);

105:   /* Allocate vectors */
106:   VecCreateSeq(PETSC_COMM_SELF, user.ndim, &user.y);
107:   VecDuplicate(user.y, &user.s);
108:   VecDuplicate(user.y, &x);

110:   /* Create TAO solver and set desired solution method */
111:   TaoCreate(PETSC_COMM_SELF, &tao);
112:   TaoSetType(tao, TAOLMVM);

114:   /* Set solution vector with an initial guess */
115:   FormInitialGuess(&user, x);
116:   TaoSetSolution(tao, x);

118:   /* Set routine for function and gradient evaluation */
119:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);

121:   /* From command line options, determine if using matrix-free hessian */
122:   PetscOptionsHasName(NULL, NULL, "-my_tao_mf", &flg);
123:   if (flg) {
124:     MatCreateShell(PETSC_COMM_SELF, user.ndim, user.ndim, user.ndim, user.ndim, (void *)&user, &H);
125:     MatShellSetOperation(H, MATOP_MULT, (void (*)(void))HessianProductMat);
126:     MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE);

128:     TaoSetHessian(tao, H, H, MatrixFreeHessian, (void *)&user);
129:   } else {
130:     MatCreateSeqAIJ(PETSC_COMM_SELF, user.ndim, user.ndim, 5, NULL, &H);
131:     MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE);
132:     TaoSetHessian(tao, H, H, FormHessian, (void *)&user);
133:   }

135:   /* Test the LMVM matrix */
136:   if (test_lmvm) {
137:     PetscOptionsSetValue(NULL, "-tao_type", "bntr");
138:     PetscOptionsSetValue(NULL, "-tao_bnk_pc_type", "lmvm");
139:   }

141:   /* Check for any TAO command line options */
142:   TaoSetFromOptions(tao);

144:   /* SOLVE THE APPLICATION */
145:   TaoSolve(tao);

147:   /* Test the LMVM matrix */
148:   if (test_lmvm) {
149:     TaoGetKSP(tao, &ksp);
150:     KSPGetPC(ksp, &pc);
151:     PCLMVMGetMatLMVM(pc, &M);
152:     VecDuplicate(x, &in);
153:     VecDuplicate(x, &out);
154:     VecDuplicate(x, &out2);
155:     VecSet(in, 5.0);
156:     MatMult(M, in, out);
157:     MatSolve(M, out, out2);
158:     VecAXPY(out2, -1.0, in);
159:     VecNorm(out2, NORM_2, &mult_solve_dist);
160:     PetscPrintf(PetscObjectComm((PetscObject)tao), "error between MatMult and MatSolve: %e\n", (double)mult_solve_dist);
161:     VecDestroy(&in);
162:     VecDestroy(&out);
163:     VecDestroy(&out2);
164:   }

166:   TaoDestroy(&tao);
167:   VecDestroy(&user.s);
168:   VecDestroy(&user.y);
169:   VecDestroy(&x);
170:   MatDestroy(&H);

172:   PetscFinalize();
173:   return 0;
174: }

176: /* ------------------------------------------------------------------- */
177: /*
178:     FormInitialGuess - Computes an initial approximation to the solution.

180:     Input Parameters:
181: .   user - user-defined application context
182: .   X    - vector

184:     Output Parameters:
185: .   X    - vector
186: */
187: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
188: {
189:   PetscReal hx = user->hx, hy = user->hy, temp;
190:   PetscReal val;
191:   PetscInt  i, j, k, nx = user->mx, ny = user->my;

193:   /* Compute initial guess */
195:   for (j = 0; j < ny; j++) {
196:     temp = PetscMin(j + 1, ny - j) * hy;
197:     for (i = 0; i < nx; i++) {
198:       k   = nx * j + i;
199:       val = PetscMin((PetscMin(i + 1, nx - i)) * hx, temp);
200:       VecSetValues(X, 1, &k, &val, ADD_VALUES);
201:     }
202:   }
203:   VecAssemblyBegin(X);
204:   VecAssemblyEnd(X);
205:   return 0;
206: }

208: /* ------------------------------------------------------------------- */
209: /*
210:    FormFunctionGradient - Evaluates the function and corresponding gradient.

212:    Input Parameters:
213:    tao - the Tao context
214:    X   - the input vector
215:    ptr - optional user-defined context, as set by TaoSetFunction()

217:    Output Parameters:
218:    f   - the newly evaluated function
219:    G   - the newly evaluated gradient
220: */
221: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr)
222: {
224:   FormFunction(tao, X, f, ptr);
225:   FormGradient(tao, X, G, ptr);
226:   return 0;
227: }

229: /* ------------------------------------------------------------------- */
230: /*
231:    FormFunction - Evaluates the function, f(X).

233:    Input Parameters:
234: .  tao - the Tao context
235: .  X   - the input vector
236: .  ptr - optional user-defined context, as set by TaoSetFunction()

238:    Output Parameters:
239: .  f    - the newly evaluated function
240: */
241: PetscErrorCode FormFunction(Tao tao, Vec X, PetscReal *f, void *ptr)
242: {
243:   AppCtx            *user = (AppCtx *)ptr;
244:   PetscReal          hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
245:   PetscReal          zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
246:   PetscReal          v, cdiv3 = user->param / three;
247:   const PetscScalar *x;
248:   PetscInt           nx = user->mx, ny = user->my, i, j, k;

251:   /* Get pointer to vector data */
252:   VecGetArrayRead(X, &x);

254:   /* Compute function contributions over the lower triangular elements */
255:   for (j = -1; j < ny; j++) {
256:     for (i = -1; i < nx; i++) {
257:       k  = nx * j + i;
258:       v  = zero;
259:       vr = zero;
260:       vt = zero;
261:       if (i >= 0 && j >= 0) v = x[k];
262:       if (i < nx - 1 && j > -1) vr = x[k + 1];
263:       if (i > -1 && j < ny - 1) vt = x[k + nx];
264:       dvdx = (vr - v) / hx;
265:       dvdy = (vt - v) / hy;
266:       fquad += dvdx * dvdx + dvdy * dvdy;
267:       flin -= cdiv3 * (v + vr + vt);
268:     }
269:   }

271:   /* Compute function contributions over the upper triangular elements */
272:   for (j = 0; j <= ny; j++) {
273:     for (i = 0; i <= nx; i++) {
274:       k  = nx * j + i;
275:       vb = zero;
276:       vl = zero;
277:       v  = zero;
278:       if (i < nx && j > 0) vb = x[k - nx];
279:       if (i > 0 && j < ny) vl = x[k - 1];
280:       if (i < nx && j < ny) v = x[k];
281:       dvdx  = (v - vl) / hx;
282:       dvdy  = (v - vb) / hy;
283:       fquad = fquad + dvdx * dvdx + dvdy * dvdy;
284:       flin  = flin - cdiv3 * (vb + vl + v);
285:     }
286:   }

288:   /* Restore vector */
289:   VecRestoreArrayRead(X, &x);

291:   /* Scale the function */
292:   area = p5 * hx * hy;
293:   *f   = area * (p5 * fquad + flin);

295:   PetscLogFlops(24.0 * nx * ny);
296:   return 0;
297: }

299: /* ------------------------------------------------------------------- */
300: /*
301:     FormGradient - Evaluates the gradient, G(X).

303:     Input Parameters:
304: .   tao  - the Tao context
305: .   X    - input vector
306: .   ptr  - optional user-defined context

308:     Output Parameters:
309: .   G - vector containing the newly evaluated gradient
310: */
311: PetscErrorCode FormGradient(Tao tao, Vec X, Vec G, void *ptr)
312: {
313:   AppCtx            *user = (AppCtx *)ptr;
314:   PetscReal          zero = 0.0, p5 = 0.5, three = 3.0, area, val;
315:   PetscInt           nx = user->mx, ny = user->my, ind, i, j, k;
316:   PetscReal          hx = user->hx, hy = user->hy;
317:   PetscReal          vb, vl, vr, vt, dvdx, dvdy;
318:   PetscReal          v, cdiv3 = user->param / three;
319:   const PetscScalar *x;

322:   /* Initialize gradient to zero */
323:   VecSet(G, zero);

325:   /* Get pointer to vector data */
326:   VecGetArrayRead(X, &x);

328:   /* Compute gradient contributions over the lower triangular elements */
329:   for (j = -1; j < ny; j++) {
330:     for (i = -1; i < nx; i++) {
331:       k  = nx * j + i;
332:       v  = zero;
333:       vr = zero;
334:       vt = zero;
335:       if (i >= 0 && j >= 0) v = x[k];
336:       if (i < nx - 1 && j > -1) vr = x[k + 1];
337:       if (i > -1 && j < ny - 1) vt = x[k + nx];
338:       dvdx = (vr - v) / hx;
339:       dvdy = (vt - v) / hy;
340:       if (i != -1 && j != -1) {
341:         ind = k;
342:         val = -dvdx / hx - dvdy / hy - cdiv3;
343:         VecSetValues(G, 1, &ind, &val, ADD_VALUES);
344:       }
345:       if (i != nx - 1 && j != -1) {
346:         ind = k + 1;
347:         val = dvdx / hx - cdiv3;
348:         VecSetValues(G, 1, &ind, &val, ADD_VALUES);
349:       }
350:       if (i != -1 && j != ny - 1) {
351:         ind = k + nx;
352:         val = dvdy / hy - cdiv3;
353:         VecSetValues(G, 1, &ind, &val, ADD_VALUES);
354:       }
355:     }
356:   }

358:   /* Compute gradient contributions over the upper triangular elements */
359:   for (j = 0; j <= ny; j++) {
360:     for (i = 0; i <= nx; i++) {
361:       k  = nx * j + i;
362:       vb = zero;
363:       vl = zero;
364:       v  = zero;
365:       if (i < nx && j > 0) vb = x[k - nx];
366:       if (i > 0 && j < ny) vl = x[k - 1];
367:       if (i < nx && j < ny) v = x[k];
368:       dvdx = (v - vl) / hx;
369:       dvdy = (v - vb) / hy;
370:       if (i != nx && j != 0) {
371:         ind = k - nx;
372:         val = -dvdy / hy - cdiv3;
373:         VecSetValues(G, 1, &ind, &val, ADD_VALUES);
374:       }
375:       if (i != 0 && j != ny) {
376:         ind = k - 1;
377:         val = -dvdx / hx - cdiv3;
378:         VecSetValues(G, 1, &ind, &val, ADD_VALUES);
379:       }
380:       if (i != nx && j != ny) {
381:         ind = k;
382:         val = dvdx / hx + dvdy / hy - cdiv3;
383:         VecSetValues(G, 1, &ind, &val, ADD_VALUES);
384:       }
385:     }
386:   }
387:   VecRestoreArrayRead(X, &x);

389:   /* Assemble gradient vector */
390:   VecAssemblyBegin(G);
391:   VecAssemblyEnd(G);

393:   /* Scale the gradient */
394:   area = p5 * hx * hy;
395:   VecScale(G, area);
396:   PetscLogFlops(24.0 * nx * ny);
397:   return 0;
398: }

400: /* ------------------------------------------------------------------- */
401: /*
402:    FormHessian - Forms the Hessian matrix.

404:    Input Parameters:
405: .  tao - the Tao context
406: .  X    - the input vector
407: .  ptr  - optional user-defined context, as set by TaoSetHessian()

409:    Output Parameters:
410: .  H     - Hessian matrix
411: .  PrecH - optionally different preconditioning Hessian
412: .  flag  - flag indicating matrix structure

414:    Notes:
415:    This routine is intended simply as an example of the interface
416:    to a Hessian evaluation routine.  Since this example compute the
417:    Hessian a column at a time, it is not particularly efficient and
418:    is not recommended.
419: */
420: PetscErrorCode FormHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
421: {
422:   AppCtx    *user = (AppCtx *)ptr;
423:   PetscInt   i, j, ndim = user->ndim;
424:   PetscReal *y, zero = 0.0, one = 1.0;
425:   PetscBool  assembled;

428:   user->xvec = X;

430:   /* Initialize Hessian entries and work vector to zero */
431:   MatAssembled(H, &assembled);
432:   if (assembled) MatZeroEntries(H);

434:   VecSet(user->s, zero);

436:   /* Loop over matrix columns to compute entries of the Hessian */
437:   for (j = 0; j < ndim; j++) {
438:     VecSetValues(user->s, 1, &j, &one, INSERT_VALUES);
439:     VecAssemblyBegin(user->s);
440:     VecAssemblyEnd(user->s);

442:     HessianProduct(ptr, user->s, user->y);

444:     VecSetValues(user->s, 1, &j, &zero, INSERT_VALUES);
445:     VecAssemblyBegin(user->s);
446:     VecAssemblyEnd(user->s);

448:     VecGetArray(user->y, &y);
449:     for (i = 0; i < ndim; i++) {
450:       if (y[i] != zero) MatSetValues(H, 1, &i, 1, &j, &y[i], ADD_VALUES);
451:     }
452:     VecRestoreArray(user->y, &y);
453:   }
454:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
455:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
456:   return 0;
457: }

459: /* ------------------------------------------------------------------- */
460: /*
461:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
462:    products.

464:    Input Parameters:
465: .  tao - the Tao context
466: .  X    - the input vector
467: .  ptr  - optional user-defined context, as set by TaoSetHessian()

469:    Output Parameters:
470: .  H     - Hessian matrix
471: .  PrecH - optionally different preconditioning Hessian
472: .  flag  - flag indicating matrix structure
473: */
474: PetscErrorCode MatrixFreeHessian(Tao tao, Vec X, Mat H, Mat PrecH, void *ptr)
475: {
476:   AppCtx *user = (AppCtx *)ptr;

478:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
480:   user->xvec = X;
481:   return 0;
482: }

484: /* ------------------------------------------------------------------- */
485: /*
486:    HessianProductMat - Computes the matrix-vector product
487:    y = mat*svec.

489:    Input Parameters:
490: .  mat  - input matrix
491: .  svec - input vector

493:    Output Parameters:
494: .  y    - solution vector
495: */
496: PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
497: {
498:   void *ptr;

501:   MatShellGetContext(mat, &ptr);
502:   HessianProduct(ptr, svec, y);
503:   return 0;
504: }

506: /* ------------------------------------------------------------------- */
507: /*
508:    Hessian Product - Computes the matrix-vector product:
509:    y = f''(x)*svec.

511:    Input Parameters:
512: .  ptr  - pointer to the user-defined context
513: .  svec - input vector

515:    Output Parameters:
516: .  y    - product vector
517: */
518: PetscErrorCode HessianProduct(void *ptr, Vec svec, Vec y)
519: {
520:   AppCtx            *user = (AppCtx *)ptr;
521:   PetscReal          p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
522:   const PetscScalar *x, *s;
523:   PetscReal          v, vb, vl, vr, vt, hxhx, hyhy;
524:   PetscInt           nx, ny, i, j, k, ind;

527:   nx   = user->mx;
528:   ny   = user->my;
529:   hx   = user->hx;
530:   hy   = user->hy;
531:   hxhx = one / (hx * hx);
532:   hyhy = one / (hy * hy);

534:   /* Get pointers to vector data */
535:   VecGetArrayRead(user->xvec, &x);
536:   VecGetArrayRead(svec, &s);

538:   /* Initialize product vector to zero */
539:   VecSet(y, zero);

541:   /* Compute f''(x)*s over the lower triangular elements */
542:   for (j = -1; j < ny; j++) {
543:     for (i = -1; i < nx; i++) {
544:       k  = nx * j + i;
545:       v  = zero;
546:       vr = zero;
547:       vt = zero;
548:       if (i != -1 && j != -1) v = s[k];
549:       if (i != nx - 1 && j != -1) {
550:         vr  = s[k + 1];
551:         ind = k + 1;
552:         val = hxhx * (vr - v);
553:         VecSetValues(y, 1, &ind, &val, ADD_VALUES);
554:       }
555:       if (i != -1 && j != ny - 1) {
556:         vt  = s[k + nx];
557:         ind = k + nx;
558:         val = hyhy * (vt - v);
559:         VecSetValues(y, 1, &ind, &val, ADD_VALUES);
560:       }
561:       if (i != -1 && j != -1) {
562:         ind = k;
563:         val = hxhx * (v - vr) + hyhy * (v - vt);
564:         VecSetValues(y, 1, &ind, &val, ADD_VALUES);
565:       }
566:     }
567:   }

569:   /* Compute f''(x)*s over the upper triangular elements */
570:   for (j = 0; j <= ny; j++) {
571:     for (i = 0; i <= nx; i++) {
572:       k  = nx * j + i;
573:       v  = zero;
574:       vl = zero;
575:       vb = zero;
576:       if (i != nx && j != ny) v = s[k];
577:       if (i != nx && j != 0) {
578:         vb  = s[k - nx];
579:         ind = k - nx;
580:         val = hyhy * (vb - v);
581:         VecSetValues(y, 1, &ind, &val, ADD_VALUES);
582:       }
583:       if (i != 0 && j != ny) {
584:         vl  = s[k - 1];
585:         ind = k - 1;
586:         val = hxhx * (vl - v);
587:         VecSetValues(y, 1, &ind, &val, ADD_VALUES);
588:       }
589:       if (i != nx && j != ny) {
590:         ind = k;
591:         val = hxhx * (v - vl) + hyhy * (v - vb);
592:         VecSetValues(y, 1, &ind, &val, ADD_VALUES);
593:       }
594:     }
595:   }
596:   /* Restore vector data */
597:   VecRestoreArrayRead(svec, &s);
598:   VecRestoreArrayRead(user->xvec, &x);

600:   /* Assemble vector */
601:   VecAssemblyBegin(y);
602:   VecAssemblyEnd(y);

604:   /* Scale resulting vector by area */
605:   area = p5 * hx * hy;
606:   VecScale(y, area);
607:   PetscLogFlops(18.0 * nx * ny);
608:   return 0;
609: }

611: /*TEST

613:    build:
614:       requires: !complex

616:    test:
617:       suffix: 1
618:       args: -tao_smonitor -tao_type ntl -tao_gatol 1.e-4

620:    test:
621:       suffix: 2
622:       args: -tao_smonitor -tao_type ntr -tao_gatol 1.e-4

624:    test:
625:       suffix: 3
626:       args: -tao_smonitor -tao_type bntr -tao_gatol 1.e-4 -my_tao_mf -tao_test_hessian

628:    test:
629:      suffix: 4
630:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnls

632:    test:
633:      suffix: 5
634:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type blmvm

636:    test:
637:      suffix: 6
638:      args: -tao_smonitor -tao_gatol 1e-3 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1

640: TEST*/