Actual source code: eptorsion1.c

petsc-master 2017-11-20
Report Typos and Errors
  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    - sysem 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>


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

 42: /*T
 43:    Concepts: TAO^Solving an unconstrained minimization problem
 44:    Routines: TaoCreate(); TaoSetType();
 45:    Routines: TaoSetInitialVector();
 46:    Routines: TaoSetObjectiveAndGradientRoutine();
 47:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 48:    Routines: TaoGetKSP(); TaoSolve();
 49:    Routines: TaoDestroy();
 50:    Processors: 1
 51: T*/

 53: /*
 54:    User-defined application context - contains data needed by the
 55:    application-provided call-back routines, FormFunction(),
 56:    FormGradient(), and FormHessian().
 57: */

 59: typedef struct {
 60:    PetscReal  param;      /* nonlinearity parameter */
 61:    PetscInt   mx, my;     /* discretization in x- and y-directions */
 62:    PetscInt   ndim;       /* problem dimension */
 63:    Vec        s, y, xvec; /* work space for computing Hessian */
 64:    PetscReal  hx, hy;     /* mesh spacing in x- and y-directions */
 65: } AppCtx;

 67: /* -------- User-defined Routines --------- */

 69: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
 70: PetscErrorCode FormFunction(Tao,Vec,PetscReal*,void*);
 71: PetscErrorCode FormGradient(Tao,Vec,Vec,void*);
 72: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat, void*);
 73: PetscErrorCode HessianProductMat(Mat,Vec,Vec);
 74: PetscErrorCode HessianProduct(void*,Vec,Vec);
 75: PetscErrorCode MatrixFreeHessian(Tao,Vec,Mat,Mat,void*);
 76: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal *,Vec,void *);

 78: PetscErrorCode main(int argc,char **argv)
 79: {
 80:   PetscErrorCode     ierr;                /* used to check for functions returning nonzeros */
 81:   PetscInt           mx=10;               /* discretization in x-direction */
 82:   PetscInt           my=10;               /* discretization in y-direction */
 83:   Vec                x;                   /* solution, gradient vectors */
 84:   PetscBool          flg;                 /* A return value when checking for use options */
 85:   Tao                tao;                 /* Tao solver context */
 86:   Mat                H;                   /* Hessian matrix */
 87:   AppCtx             user;                /* application context */
 88:   PetscMPIInt        size;                /* number of processes */
 89:   PetscReal          one=1.0;

 91:   /* Initialize TAO,PETSc */
 92:   PetscInitialize(&argc,&argv,(char *)0,help);
 93:   MPI_Comm_size(MPI_COMM_WORLD,&size);
 94:   if (size >1) SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");

 96:   /* Specify default parameters for the problem, check for command-line overrides */
 97:   user.param = 5.0;
 98:   PetscOptionsGetInt(NULL,NULL,"-my",&my,&flg);
 99:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,&flg);
100:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);

102:   PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
103:   PetscPrintf(PETSC_COMM_SELF,"mx: %D     my: %D   \n\n",mx,my);
104:   user.ndim = mx * my; user.mx = mx; user.my = my;
105:   user.hx = one/(mx+1); user.hy = one/(my+1);

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

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

116:   /* Set solution vector with an initial guess */
117:   FormInitialGuess(&user,x);
118:   TaoSetInitialVector(tao,x);

120:   /* Set routine for function and gradient evaluation */
121:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

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

130:     TaoSetHessianRoutine(tao,H,H,MatrixFreeHessian,(void *)&user);

132:     /* Set null preconditioner.  Alternatively, set user-provided
133:        preconditioner or explicitly form preconditioning matrix */
134:     PetscOptionsSetValue(NULL,"-pc_type","none");
135:   } else {
136:     MatCreateSeqAIJ(PETSC_COMM_SELF,user.ndim,user.ndim,5,NULL,&H);
137:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);
138:     TaoSetHessianRoutine(tao,H,H,FormHessian,(void *)&user);
139:   }


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

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

148:   TaoDestroy(&tao);
149:   VecDestroy(&user.s);
150:   VecDestroy(&user.y);
151:   VecDestroy(&x);
152:   MatDestroy(&H);

154:   PetscFinalize();
155:   return ierr;
156: }

158: /* ------------------------------------------------------------------- */
159: /*
160:     FormInitialGuess - Computes an initial approximation to the solution.

162:     Input Parameters:
163: .   user - user-defined application context
164: .   X    - vector

166:     Output Parameters:
167: .   X    - vector
168: */
169: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
170: {
171:   PetscReal      hx = user->hx, hy = user->hy, temp;
172:   PetscReal      val;
174:   PetscInt       i, j, k, nx = user->mx, ny = user->my;

176:   /* Compute initial guess */
178:   for (j=0; j<ny; j++) {
179:     temp = PetscMin(j+1,ny-j)*hy;
180:     for (i=0; i<nx; i++) {
181:       k   = nx*j + i;
182:       val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
183:       VecSetValues(X,1,&k,&val,ADD_VALUES);
184:     }
185:   }
186:   VecAssemblyBegin(X);
187:   VecAssemblyEnd(X);
188:   return(0);
189: }

191: /* ------------------------------------------------------------------- */
192: /*
193:    FormFunctionGradient - Evaluates the function and corresponding gradient.

195:    Input Parameters:
196:    tao - the Tao context
197:    X   - the input vector
198:    ptr - optional user-defined context, as set by TaoSetFunction()

200:    Output Parameters:
201:    f   - the newly evaluated function
202:    G   - the newly evaluated gradient
203: */
204: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
205: {

209:   FormFunction(tao,X,f,ptr);
210:   FormGradient(tao,X,G,ptr);
211:   return(0);
212: }

214: /* ------------------------------------------------------------------- */
215: /*
216:    FormFunction - Evaluates the function, f(X).

218:    Input Parameters:
219: .  tao - the Tao context
220: .  X   - the input vector
221: .  ptr - optional user-defined context, as set by TaoSetFunction()

223:    Output Parameters:
224: .  f    - the newly evaluated function
225: */
226: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
227: {
228:   AppCtx            *user = (AppCtx *) ptr;
229:   PetscReal         hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
230:   PetscReal         zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
231:   PetscReal         v, cdiv3 = user->param/three;
232:   const PetscScalar *x;
233:   PetscErrorCode    ierr;
234:   PetscInt          nx = user->mx, ny = user->my, i, j, k;

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

240:   /* Compute function contributions over the lower triangular elements */
241:   for (j=-1; j<ny; j++) {
242:     for (i=-1; i<nx; i++) {
243:       k = nx*j + i;
244:       v = zero;
245:       vr = zero;
246:       vt = zero;
247:       if (i >= 0 && j >= 0) v = x[k];
248:       if (i < nx-1 && j > -1) vr = x[k+1];
249:       if (i > -1 && j < ny-1) vt = x[k+nx];
250:       dvdx = (vr-v)/hx;
251:       dvdy = (vt-v)/hy;
252:       fquad += dvdx*dvdx + dvdy*dvdy;
253:       flin -= cdiv3*(v+vr+vt);
254:     }
255:   }

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

274:   /* Restore vector */
275:   VecRestoreArrayRead(X,&x);

277:   /* Scale the function */
278:   area = p5*hx*hy;
279:   *f = area*(p5*fquad+flin);

281:   PetscLogFlops(nx*ny*24);
282:   return(0);
283: }

285: /* ------------------------------------------------------------------- */
286: /*
287:     FormGradient - Evaluates the gradient, G(X).

289:     Input Parameters:
290: .   tao  - the Tao context
291: .   X    - input vector
292: .   ptr  - optional user-defined context

294:     Output Parameters:
295: .   G - vector containing the newly evaluated gradient
296: */
297: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
298: {
299:   AppCtx            *user = (AppCtx *) ptr;
300:   PetscReal         zero=0.0, p5=0.5, three = 3.0, area, val;
301:   PetscErrorCode    ierr;
302:   PetscInt          nx = user->mx, ny = user->my, ind, i, j, k;
303:   PetscReal         hx = user->hx, hy = user->hy;
304:   PetscReal         vb, vl, vr, vt, dvdx, dvdy;
305:   PetscReal         v, cdiv3 = user->param/three;
306:   const PetscScalar *x;

309:   /* Initialize gradient to zero */
310:   VecSet(G, zero);

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

315:   /* Compute gradient contributions over the lower triangular elements */
316:   for (j=-1; j<ny; j++) {
317:     for (i=-1; i<nx; i++) {
318:       k  = nx*j + i;
319:       v  = zero;
320:       vr = zero;
321:       vt = zero;
322:       if (i >= 0 && j >= 0)    v = x[k];
323:       if (i < nx-1 && j > -1) vr = x[k+1];
324:       if (i > -1 && j < ny-1) vt = x[k+nx];
325:       dvdx = (vr-v)/hx;
326:       dvdy = (vt-v)/hy;
327:       if (i != -1 && j != -1) {
328:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
329:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
330:       }
331:       if (i != nx-1 && j != -1) {
332:         ind = k+1; val =  dvdx/hx - cdiv3;
333:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
334:       }
335:       if (i != -1 && j != ny-1) {
336:         ind = k+nx; val = dvdy/hy - cdiv3;
337:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
338:       }
339:     }
340:   }

342:   /* Compute gradient contributions over the upper triangular elements */
343:   for (j=0; j<=ny; j++) {
344:     for (i=0; i<=nx; i++) {
345:       k = nx*j + i;
346:       vb = zero;
347:       vl = zero;
348:       v  = zero;
349:       if (i < nx && j > 0) vb = x[k-nx];
350:       if (i > 0 && j < ny) vl = x[k-1];
351:       if (i < nx && j < ny) v = x[k];
352:       dvdx = (v-vl)/hx;
353:       dvdy = (v-vb)/hy;
354:       if (i != nx && j != 0) {
355:         ind = k-nx; val = - dvdy/hy - cdiv3;
356:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
357:       }
358:       if (i != 0 && j != ny) {
359:         ind = k-1; val =  - dvdx/hx - cdiv3;
360:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
361:       }
362:       if (i != nx && j != ny) {
363:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
364:         VecSetValues(G,1,&ind,&val,ADD_VALUES);
365:       }
366:     }
367:   }
368:   VecRestoreArrayRead(X,&x);

370:   /* Assemble gradient vector */
371:   VecAssemblyBegin(G);
372:   VecAssemblyEnd(G);

374:   /* Scale the gradient */
375:   area = p5*hx*hy;
376:   VecScale(G, area);
377:   PetscLogFlops(nx*ny*24);
378:   return(0);
379: }

381: /* ------------------------------------------------------------------- */
382: /*
383:    FormHessian - Forms the Hessian matrix.

385:    Input Parameters:
386: .  tao - the Tao context
387: .  X    - the input vector
388: .  ptr  - optional user-defined context, as set by TaoSetHessian()

390:    Output Parameters:
391: .  H     - Hessian matrix
392: .  PrecH - optionally different preconditioning Hessian
393: .  flag  - flag indicating matrix structure

395:    Notes:
396:    This routine is intended simply as an example of the interface
397:    to a Hessian evaluation routine.  Since this example compute the
398:    Hessian a column at a time, it is not particularly efficient and
399:    is not recommended.
400: */
401: PetscErrorCode FormHessian(Tao tao,Vec X,Mat H,Mat Hpre, void *ptr)
402: {
403:   AppCtx         *user = (AppCtx *) ptr;
405:   PetscInt       i,j, ndim = user->ndim;
406:   PetscReal      *y, zero = 0.0, one = 1.0;
407:   PetscBool      assembled;

410:   user->xvec = X;

412:   /* Initialize Hessian entries and work vector to zero */
413:   MatAssembled(H,&assembled);
414:   if (assembled){MatZeroEntries(H); }

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

418:   /* Loop over matrix columns to compute entries of the Hessian */
419:   for (j=0; j<ndim; j++) {
420:     VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
421:     VecAssemblyBegin(user->s);
422:     VecAssemblyEnd(user->s);

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

426:     VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
427:     VecAssemblyBegin(user->s);
428:     VecAssemblyEnd(user->s);

430:     VecGetArray(user->y,&y);
431:     for (i=0; i<ndim; i++) {
432:       if (y[i] != zero) {
433:         MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
434:       }
435:     }
436:     VecRestoreArray(user->y,&y);
437:   }
438:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
439:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
440:   return(0);
441: }

443: /* ------------------------------------------------------------------- */
444: /*
445:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
446:    products.

448:    Input Parameters:
449: .  tao - the Tao context
450: .  X    - the input vector
451: .  ptr  - optional user-defined context, as set by TaoSetHessian()

453:    Output Parameters:
454: .  H     - Hessian matrix
455: .  PrecH - optionally different preconditioning Hessian
456: .  flag  - flag indicating matrix structure
457: */
458: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
459: {
460:   AppCtx     *user = (AppCtx *) ptr;

462:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
463:   user->xvec = X;
464:   return(0);
465: }

467: /* ------------------------------------------------------------------- */
468: /*
469:    HessianProductMat - Computes the matrix-vector product
470:    y = mat*svec.

472:    Input Parameters:
473: .  mat  - input matrix
474: .  svec - input vector

476:    Output Parameters:
477: .  y    - solution vector
478: */
479: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
480: {
481:   void           *ptr;

485:   MatShellGetContext(mat,&ptr);
486:   HessianProduct(ptr,svec,y);
487:   return(0);
488: }

490: /* ------------------------------------------------------------------- */
491: /*
492:    Hessian Product - Computes the matrix-vector product:
493:    y = f''(x)*svec.

495:    Input Parameters
496: .  ptr  - pointer to the user-defined context
497: .  svec - input vector

499:    Output Parameters:
500: .  y    - product vector
501: */
502: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
503: {
504:   AppCtx            *user = (AppCtx *)ptr;
505:   PetscReal         p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area;
506:   const PetscScalar *x, *s;
507:   PetscReal         v, vb, vl, vr, vt, hxhx, hyhy;
508:   PetscErrorCode    ierr;
509:   PetscInt          nx, ny, i, j, k, ind;

512:   nx   = user->mx;
513:   ny   = user->my;
514:   hx   = user->hx;
515:   hy   = user->hy;
516:   hxhx = one/(hx*hx);
517:   hyhy = one/(hy*hy);

519:   /* Get pointers to vector data */
520:   VecGetArrayRead(user->xvec,&x);
521:   VecGetArrayRead(svec,&s);

523:   /* Initialize product vector to zero */
524:   VecSet(y, zero);

526:   /* Compute f''(x)*s over the lower triangular elements */
527:   for (j=-1; j<ny; j++) {
528:     for (i=-1; i<nx; i++) {
529:        k = nx*j + i;
530:        v = zero;
531:        vr = zero;
532:        vt = zero;
533:        if (i != -1 && j != -1) v = s[k];
534:        if (i != nx-1 && j != -1) {
535:          vr = s[k+1];
536:          ind = k+1; val = hxhx*(vr-v);
537:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
538:        }
539:        if (i != -1 && j != ny-1) {
540:          vt = s[k+nx];
541:          ind = k+nx; val = hyhy*(vt-v);
542:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
543:        }
544:        if (i != -1 && j != -1) {
545:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
546:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
547:        }
548:      }
549:    }

551:   /* Compute f''(x)*s over the upper triangular elements */
552:   for (j=0; j<=ny; j++) {
553:     for (i=0; i<=nx; i++) {
554:        k = nx*j + i;
555:        v = zero;
556:        vl = zero;
557:        vb = zero;
558:        if (i != nx && j != ny) v = s[k];
559:        if (i != nx && j != 0) {
560:          vb = s[k-nx];
561:          ind = k-nx; val = hyhy*(vb-v);
562:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
563:        }
564:        if (i != 0 && j != ny) {
565:          vl = s[k-1];
566:          ind = k-1; val = hxhx*(vl-v);
567:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
568:        }
569:        if (i != nx && j != ny) {
570:          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
571:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
572:        }
573:     }
574:   }
575:   /* Restore vector data */
576:   VecRestoreArrayRead(svec,&s);
577:   VecRestoreArrayRead(user->xvec,&x);

579:   /* Assemble vector */
580:   VecAssemblyBegin(y);
581:   VecAssemblyEnd(y);

583:   /* Scale resulting vector by area */
584:   area = p5*hx*hy;
585:   VecScale(y, area);
586:   PetscLogFlops(nx*ny*18);
587:   return(0);
588: }