Actual source code: eptorsion1.c

petsc-3.7.3 2016-08-01
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 *);

 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:   AppCtx             user;                /* application context */
 90:   PetscMPIInt        size;                /* number of processes */
 91:   PetscReal          one=1.0;

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

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

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

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

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

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

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

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

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

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

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


145:   /* Check for any TAO command line options */
146:   TaoSetFromOptions(tao);

148:   /* SOLVE THE APPLICATION */
149:   TaoSolve(tao);

151:   TaoDestroy(&tao);
152:   VecDestroy(&user.s);
153:   VecDestroy(&user.y);
154:   VecDestroy(&x);
155:   MatDestroy(&H);

157:   PetscFinalize();
158:   return 0;
159: }

161: /* ------------------------------------------------------------------- */
164: /*
165:     FormInitialGuess - Computes an initial approximation to the solution.

167:     Input Parameters:
168: .   user - user-defined application context
169: .   X    - vector

171:     Output Parameters:
172: .   X    - vector
173: */
174: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
175: {
176:   PetscReal      hx = user->hx, hy = user->hy, temp;
177:   PetscReal      val;
179:   PetscInt       i, j, k, nx = user->mx, ny = user->my;

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

195: /* ------------------------------------------------------------------- */
198: /*
199:    FormFunctionGradient - Evaluates the function and corresponding gradient.

201:    Input Parameters:
202:    tao - the Tao context
203:    X   - the input vector
204:    ptr - optional user-defined context, as set by TaoSetFunction()

206:    Output Parameters:
207:    f   - the newly evaluated function
208:    G   - the newly evaluated gradient
209: */
210: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
211: {

214:   FormFunction(tao,X,f,ptr);
215:   FormGradient(tao,X,G,ptr);
216:   return 0;
217: }

219: /* ------------------------------------------------------------------- */
222: /*
223:    FormFunction - Evaluates the function, f(X).

225:    Input Parameters:
226: .  tao - the Tao context
227: .  X   - the input vector
228: .  ptr - optional user-defined context, as set by TaoSetFunction()

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

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

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

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

280:   /* Restore vector */
281:   VecRestoreArray(X,&x);

283:   /* Scale the function */
284:   area = p5*hx*hy;
285:   *f = area*(p5*fquad+flin);

287:   PetscLogFlops(nx*ny*24);
288:   return 0;
289: }

291: /* ------------------------------------------------------------------- */
294: /*
295:     FormGradient - Evaluates the gradient, G(X).

297:     Input Parameters:
298: .   tao  - the Tao context
299: .   X    - input vector
300: .   ptr  - optional user-defined context

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

315:   /* Initialize gradient to zero */
316:   VecSet(G, zero);

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

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

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

376:   /* Assemble gradient vector */
377:   VecAssemblyBegin(G);
378:   VecAssemblyEnd(G);

380:   /* Scale the gradient */
381:   area = p5*hx*hy;
382:   VecScale(G, area);
383:   PetscLogFlops(nx*ny*24);
384:   return 0;
385: }

387: /* ------------------------------------------------------------------- */
390: /*
391:    FormHessian - Forms the Hessian matrix.

393:    Input Parameters:
394: .  tao - the Tao context
395: .  X    - the input vector
396: .  ptr  - optional user-defined context, as set by TaoSetHessian()

398:    Output Parameters:
399: .  H     - Hessian matrix
400: .  PrecH - optionally different preconditioning Hessian
401: .  flag  - flag indicating matrix structure

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

417:   user->xvec = X;

419:   /* Initialize Hessian entries and work vector to zero */
420:   MatAssembled(H,&assembled);
421:   if (assembled){MatZeroEntries(H); }

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

425:   /* Loop over matrix columns to compute entries of the Hessian */
426:   for (j=0; j<ndim; j++) {
427:     VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
428:     VecAssemblyBegin(user->s);
429:     VecAssemblyEnd(user->s);

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

433:     VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
434:     VecAssemblyBegin(user->s);
435:     VecAssemblyEnd(user->s);

437:     VecGetArray(user->y,&y);
438:     for (i=0; i<ndim; i++) {
439:       if (y[i] != zero) {
440:         MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
441:       }
442:     }
443:     VecRestoreArray(user->y,&y);
444:   }
445:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
446:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
447:   return 0;
448: }

450: /* ------------------------------------------------------------------- */
453: /*
454:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
455:    products.

457:    Input Parameters:
458: .  tao - the Tao context
459: .  X    - the input vector
460: .  ptr  - optional user-defined context, as set by TaoSetHessian()

462:    Output Parameters:
463: .  H     - Hessian matrix
464: .  PrecH - optionally different preconditioning Hessian
465: .  flag  - flag indicating matrix structure
466: */
467: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
468: {
469:   AppCtx     *user = (AppCtx *) ptr;

471:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
472:   user->xvec = X;
473:   return 0;
474: }

476: /* ------------------------------------------------------------------- */
479: /*
480:    HessianProductMat - Computes the matrix-vector product
481:    y = mat*svec.

483:    Input Parameters:
484: .  mat  - input matrix
485: .  svec - input vector

487:    Output Parameters:
488: .  y    - solution vector
489: */
490: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
491: {
492:   void           *ptr;

495:   MatShellGetContext(mat,&ptr);
496:   HessianProduct(ptr,svec,y);
497:   return 0;
498: }

500: /* ------------------------------------------------------------------- */
503: /*
504:    Hessian Product - Computes the matrix-vector product:
505:    y = f''(x)*svec.

507:    Input Parameters
508: .  ptr  - pointer to the user-defined context
509: .  svec - input vector

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

522:   nx   = user->mx;
523:   ny   = user->my;
524:   hx   = user->hx;
525:   hy   = user->hy;
526:   hxhx = one/(hx*hx);
527:   hyhy = one/(hy*hy);

529:   /* Get pointers to vector data */
530:   VecGetArray(user->xvec,&x);
531:   VecGetArray(svec,&s);

533:   /* Initialize product vector to zero */
534:   VecSet(y, zero);

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

561:   /* Compute f''(x)*s over the upper triangular elements */
562:   for (j=0; j<=ny; j++) {
563:     for (i=0; i<=nx; i++) {
564:        k = nx*j + i;
565:        v = zero;
566:        vl = zero;
567:        vb = zero;
568:        if (i != nx && j != ny) v = s[k];
569:        if (i != nx && j != 0) {
570:          vb = s[k-nx];
571:          ind = k-nx; val = hyhy*(vb-v);
572:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
573:        }
574:        if (i != 0 && j != ny) {
575:          vl = s[k-1];
576:          ind = k-1; val = hxhx*(vl-v);
577:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
578:        }
579:        if (i != nx && j != ny) {
580:          ind = k; val = hxhx*(v-vl) + hyhy*(v-vb);
581:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
582:        }
583:     }
584:   }
585:   /* Restore vector data */
586:   VecRestoreArray(svec,&s);
587:   VecRestoreArray(user->xvec,&x);

589:   /* Assemble vector */
590:   VecAssemblyBegin(y);
591:   VecAssemblyEnd(y);

593:   /* Scale resulting vector by area */
594:   area = p5*hx*hy;
595:   VecScale(y, area);
596:   PetscLogFlops(nx*ny*18);
597:   return 0;
598: }