Actual source code: eptorsion1.c

petsc-dev 2014-07-25
Report Typos and Errors
  1: /* Program usage: mpirun -np 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: TaoGetConvergedReason(); 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 *);

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

 95:   /* Initialize TAO,PETSc */
 96:   PetscInitialize(&argc,&argv,(char *)0,help);

 98:   MPI_Comm_size(MPI_COMM_WORLD,&size);
 99:   if (size >1) SETERRQ(PETSC_COMM_SELF,1,"Incorrect number of processors");

101:   /* Specify default parameters for the problem, check for command-line overrides */
102:   user.param = 5.0;
103:   PetscOptionsGetInt(NULL,"-my",&my,&flg);
104:   PetscOptionsGetInt(NULL,"-mx",&mx,&flg);
105:   PetscOptionsGetReal(NULL,"-par",&user.param,&flg);

107:   PetscPrintf(PETSC_COMM_SELF,"\n---- Elastic-Plastic Torsion Problem -----\n");
108:   PetscPrintf(PETSC_COMM_SELF,"mx: %D     my: %D   \n\n",mx,my);
109:   user.ndim = mx * my; user.mx = mx; user.my = my;
110:   user.hx = one/(mx+1); user.hy = one/(my+1);

112:   /* Allocate vectors */
113:   VecCreateSeq(PETSC_COMM_SELF,user.ndim,&user.y);
114:   VecDuplicate(user.y,&user.s);
115:   VecDuplicate(user.y,&x);

117:   /* Create TAO solver and set desired solution method */
118:   TaoCreate(PETSC_COMM_SELF,&tao);
119:   TaoSetType(tao,TAOLMVM);

121:   /* Set solution vector with an initial guess */
122:   FormInitialGuess(&user,x);
123:   TaoSetInitialVector(tao,x);

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

128:   /* From command line options, determine if using matrix-free hessian */
129:   PetscOptionsHasName(NULL,"-my_tao_mf",&flg);
130:   if (flg) {
131:     MatCreateShell(PETSC_COMM_SELF,user.ndim,user.ndim,user.ndim,user.ndim,(void*)&user,&H);
132:     MatShellSetOperation(H,MATOP_MULT,(void(*)(void))HessianProductMat);
133:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);

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

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


147:   /* Check for any TAO command line options */
148:   TaoSetFromOptions(tao);

150:   /* SOLVE THE APPLICATION */
151:   TaoSolve(tao);
152:   TaoGetKSP(tao,&ksp);
153:   if (ksp) {
154:     KSPView(ksp,PETSC_VIEWER_STDOUT_SELF);
155:   }

157:   /*
158:      To View TAO solver information use
159:       TaoView(tao);
160:   */

162:   /* Get information on termination */
163:   TaoGetConvergedReason(tao,&reason);
164:   if (reason <= 0){
165:     PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
166:   }

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

174:   PetscFinalize();
175:   return 0;
176: }

178: /* ------------------------------------------------------------------- */
181: /*
182:     FormInitialGuess - Computes an initial approximation to the solution.

184:     Input Parameters:
185: .   user - user-defined application context
186: .   X    - vector

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

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

212: /* ------------------------------------------------------------------- */
215: /*
216:    FormFunctionGradient - Evaluates the function and corresponding gradient.

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:    G   - the newly evaluated gradient
226: */
227: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
228: {

231:   FormFunction(tao,X,f,ptr);
232:   FormGradient(tao,X,G,ptr);
233:   return 0;
234: }

236: /* ------------------------------------------------------------------- */
239: /*
240:    FormFunction - Evaluates the function, f(X).

242:    Input Parameters:
243: .  tao - the Tao context
244: .  X   - the input vector
245: .  ptr - optional user-defined context, as set by TaoSetFunction()

247:    Output Parameters:
248: .  f    - the newly evaluated function
249: */
250: PetscErrorCode FormFunction(Tao tao,Vec X,PetscReal *f,void *ptr)
251: {
252:   AppCtx         *user = (AppCtx *) ptr;
253:   PetscReal      hx = user->hx, hy = user->hy, area, three = 3.0, p5 = 0.5;
254:   PetscReal      zero = 0.0, vb, vl, vr, vt, dvdx, dvdy, flin = 0.0, fquad = 0.0;
255:   PetscReal      v, cdiv3 = user->param/three;
256:   PetscReal      *x;
258:   PetscInt       nx = user->mx, ny = user->my, i, j, k;

260:   /* Get pointer to vector data */
261:   VecGetArray(X,&x);

263:   /* Compute function contributions over the lower triangular elements */
264:   for (j=-1; j<ny; j++) {
265:     for (i=-1; i<nx; i++) {
266:       k = nx*j + i;
267:       v = zero;
268:       vr = zero;
269:       vt = zero;
270:       if (i >= 0 && j >= 0) v = x[k];
271:       if (i < nx-1 && j > -1) vr = x[k+1];
272:       if (i > -1 && j < ny-1) vt = x[k+nx];
273:       dvdx = (vr-v)/hx;
274:       dvdy = (vt-v)/hy;
275:       fquad += dvdx*dvdx + dvdy*dvdy;
276:       flin -= cdiv3*(v+vr+vt);
277:     }
278:   }

280:   /* Compute function contributions over the upper triangular elements */
281:   for (j=0; j<=ny; j++) {
282:     for (i=0; i<=nx; i++) {
283:       k = nx*j + i;
284:       vb = zero;
285:       vl = zero;
286:       v = zero;
287:       if (i < nx && j > 0) vb = x[k-nx];
288:       if (i > 0 && j < ny) vl = x[k-1];
289:       if (i < nx && j < ny) v = x[k];
290:       dvdx = (v-vl)/hx;
291:       dvdy = (v-vb)/hy;
292:       fquad = fquad + dvdx*dvdx + dvdy*dvdy;
293:       flin = flin - cdiv3*(vb+vl+v);
294:     }
295:   }

297:   /* Restore vector */
298:   VecRestoreArray(X,&x);

300:   /* Scale the function */
301:   area = p5*hx*hy;
302:   *f = area*(p5*fquad+flin);

304:   PetscLogFlops(nx*ny*24);
305:   return 0;
306: }

308: /* ------------------------------------------------------------------- */
311: /*
312:     FormGradient - Evaluates the gradient, G(X).

314:     Input Parameters:
315: .   tao  - the Tao context
316: .   X    - input vector
317: .   ptr  - optional user-defined context

319:     Output Parameters:
320: .   G - vector containing the newly evaluated gradient
321: */
322: PetscErrorCode FormGradient(Tao tao,Vec X,Vec G,void *ptr)
323: {
324:   AppCtx         *user = (AppCtx *) ptr;
325:   PetscReal      zero=0.0, p5=0.5, three = 3.0, area, val, *x;
327:   PetscInt       nx = user->mx, ny = user->my, ind, i, j, k;
328:   PetscReal      hx = user->hx, hy = user->hy;
329:   PetscReal      vb, vl, vr, vt, dvdx, dvdy;
330:   PetscReal      v, cdiv3 = user->param/three;

332:   /* Initialize gradient to zero */
333:   VecSet(G, zero);

335:   /* Get pointer to vector data */
336:   VecGetArray(X,&x);

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

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

393:   /* Assemble gradient vector */
394:   VecAssemblyBegin(G);
395:   VecAssemblyEnd(G);

397:   /* Scale the gradient */
398:   area = p5*hx*hy;
399:   VecScale(G, area);
400:   PetscLogFlops(nx*ny*24);
401:   return 0;
402: }

404: /* ------------------------------------------------------------------- */
407: /*
408:    FormHessian - Forms the Hessian matrix.

410:    Input Parameters:
411: .  tao - the Tao context
412: .  X    - the input vector
413: .  ptr  - optional user-defined context, as set by TaoSetHessian()

415:    Output Parameters:
416: .  H     - Hessian matrix
417: .  PrecH - optionally different preconditioning Hessian
418: .  flag  - flag indicating matrix structure

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

434:   user->xvec = X;

436:   /* Initialize Hessian entries and work vector to zero */
437:   MatAssembled(H,&assembled);
438:   if (assembled){MatZeroEntries(H); }

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

442:   /* Loop over matrix columns to compute entries of the Hessian */
443:   for (j=0; j<ndim; j++) {
444:     VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
445:     VecAssemblyBegin(user->s);
446:     VecAssemblyEnd(user->s);

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

450:     VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
451:     VecAssemblyBegin(user->s);
452:     VecAssemblyEnd(user->s);

454:     VecGetArray(user->y,&y);
455:     for (i=0; i<ndim; i++) {
456:       if (y[i] != zero) {
457:         MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
458:       }
459:     }
460:     VecRestoreArray(user->y,&y);
461:   }
462:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
463:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
464:   return 0;
465: }

467: /* ------------------------------------------------------------------- */
470: /*
471:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
472:    products.

474:    Input Parameters:
475: .  tao - the Tao context
476: .  X    - the input vector
477: .  ptr  - optional user-defined context, as set by TaoSetHessian()

479:    Output Parameters:
480: .  H     - Hessian matrix
481: .  PrecH - optionally different preconditioning Hessian
482: .  flag  - flag indicating matrix structure
483: */
484: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
485: {
486:   AppCtx     *user = (AppCtx *) ptr;

488:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
489:   user->xvec = X;
490:   return 0;
491: }

493: /* ------------------------------------------------------------------- */
496: /*
497:    HessianProductMat - Computes the matrix-vector product
498:    y = mat*svec.

500:    Input Parameters:
501: .  mat  - input matrix
502: .  svec - input vector

504:    Output Parameters:
505: .  y    - solution vector
506: */
507: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
508: {
509:   void           *ptr;

512:   MatShellGetContext(mat,&ptr);
513:   HessianProduct(ptr,svec,y);
514:   return 0;
515: }

517: /* ------------------------------------------------------------------- */
520: /*
521:    Hessian Product - Computes the matrix-vector product:
522:    y = f''(x)*svec.

524:    Input Parameters
525: .  ptr  - pointer to the user-defined context
526: .  svec - input vector

528:    Output Parameters:
529: .  y    - product vector
530: */
531: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
532: {
533:   AppCtx         *user = (AppCtx *)ptr;
534:   PetscReal      p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
535:   PetscReal      v, vb, vl, vr, vt, hxhx, hyhy;
537:   PetscInt       nx, ny, i, j, k, ind;

539:   nx   = user->mx;
540:   ny   = user->my;
541:   hx   = user->hx;
542:   hy   = user->hy;
543:   hxhx = one/(hx*hx);
544:   hyhy = one/(hy*hy);

546:   /* Get pointers to vector data */
547:   VecGetArray(user->xvec,&x);
548:   VecGetArray(svec,&s);

550:   /* Initialize product vector to zero */
551:   VecSet(y, zero);

553:   /* Compute f''(x)*s over the lower triangular elements */
554:   for (j=-1; j<ny; j++) {
555:     for (i=-1; i<nx; i++) {
556:        k = nx*j + i;
557:        v = zero;
558:        vr = zero;
559:        vt = zero;
560:        if (i != -1 && j != -1) v = s[k];
561:        if (i != nx-1 && j != -1) {
562:          vr = s[k+1];
563:          ind = k+1; val = hxhx*(vr-v);
564:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
565:        }
566:        if (i != -1 && j != ny-1) {
567:          vt = s[k+nx];
568:          ind = k+nx; val = hyhy*(vt-v);
569:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
570:        }
571:        if (i != -1 && j != -1) {
572:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
573:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
574:        }
575:      }
576:    }

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

606:   /* Assemble vector */
607:   VecAssemblyBegin(y);
608:   VecAssemblyEnd(y);

610:   /* Scale resulting vector by area */
611:   area = p5*hx*hy;
612:   VecScale(y, area);
613:   PetscLogFlops(nx*ny*18);
614:   return 0;
615: }