Actual source code: eptorsion1.c

petsc-master 2017-09-21
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 */
177:   for (j=0; j<ny; j++) {
178:     temp = PetscMin(j+1,ny-j)*hy;
179:     for (i=0; i<nx; i++) {
180:       k   = nx*j + i;
181:       val = PetscMin((PetscMin(i+1,nx-i))*hx,temp);
182:       VecSetValues(X,1,&k,&val,ADD_VALUES);
183:     }
184:   }
185:   VecAssemblyBegin(X);
186:   VecAssemblyEnd(X);
187:   return 0;
188: }

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

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

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

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

212: /* ------------------------------------------------------------------- */
213: /*
214:    FormFunction - Evaluates the function, f(X).

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

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

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

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

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

271:   /* Restore vector */
272:   VecRestoreArray(X,&x);

274:   /* Scale the function */
275:   area = p5*hx*hy;
276:   *f = area*(p5*fquad+flin);

278:   PetscLogFlops(nx*ny*24);
279:   return 0;
280: }

282: /* ------------------------------------------------------------------- */
283: /*
284:     FormGradient - Evaluates the gradient, G(X).

286:     Input Parameters:
287: .   tao  - the Tao context
288: .   X    - input vector
289: .   ptr  - optional user-defined context

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

304:   /* Initialize gradient to zero */
305:   VecSet(G, zero);

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

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

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

365:   /* Assemble gradient vector */
366:   VecAssemblyBegin(G);
367:   VecAssemblyEnd(G);

369:   /* Scale the gradient */
370:   area = p5*hx*hy;
371:   VecScale(G, area);
372:   PetscLogFlops(nx*ny*24);
373:   return 0;
374: }

376: /* ------------------------------------------------------------------- */
377: /*
378:    FormHessian - Forms the Hessian matrix.

380:    Input Parameters:
381: .  tao - the Tao context
382: .  X    - the input vector
383: .  ptr  - optional user-defined context, as set by TaoSetHessian()

385:    Output Parameters:
386: .  H     - Hessian matrix
387: .  PrecH - optionally different preconditioning Hessian
388: .  flag  - flag indicating matrix structure

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

404:   user->xvec = X;

406:   /* Initialize Hessian entries and work vector to zero */
407:   MatAssembled(H,&assembled);
408:   if (assembled){MatZeroEntries(H); }

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

412:   /* Loop over matrix columns to compute entries of the Hessian */
413:   for (j=0; j<ndim; j++) {
414:     VecSetValues(user->s,1,&j,&one,INSERT_VALUES);
415:     VecAssemblyBegin(user->s);
416:     VecAssemblyEnd(user->s);

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

420:     VecSetValues(user->s,1,&j,&zero,INSERT_VALUES);
421:     VecAssemblyBegin(user->s);
422:     VecAssemblyEnd(user->s);

424:     VecGetArray(user->y,&y);
425:     for (i=0; i<ndim; i++) {
426:       if (y[i] != zero) {
427:         MatSetValues(H,1,&i,1,&j,&y[i],ADD_VALUES);
428:       }
429:     }
430:     VecRestoreArray(user->y,&y);
431:   }
432:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
433:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
434:   return 0;
435: }

437: /* ------------------------------------------------------------------- */
438: /*
439:    MatrixFreeHessian - Sets a pointer for use in computing Hessian-vector
440:    products.

442:    Input Parameters:
443: .  tao - the Tao context
444: .  X    - the input vector
445: .  ptr  - optional user-defined context, as set by TaoSetHessian()

447:    Output Parameters:
448: .  H     - Hessian matrix
449: .  PrecH - optionally different preconditioning Hessian
450: .  flag  - flag indicating matrix structure
451: */
452: PetscErrorCode MatrixFreeHessian(Tao tao,Vec X,Mat H,Mat PrecH, void *ptr)
453: {
454:   AppCtx     *user = (AppCtx *) ptr;

456:   /* Sets location of vector for use in computing matrix-vector products  of the form H(X)*y  */
457:   user->xvec = X;
458:   return 0;
459: }

461: /* ------------------------------------------------------------------- */
462: /*
463:    HessianProductMat - Computes the matrix-vector product
464:    y = mat*svec.

466:    Input Parameters:
467: .  mat  - input matrix
468: .  svec - input vector

470:    Output Parameters:
471: .  y    - solution vector
472: */
473: PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y)
474: {
475:   void           *ptr;

478:   MatShellGetContext(mat,&ptr);
479:   HessianProduct(ptr,svec,y);
480:   return 0;
481: }

483: /* ------------------------------------------------------------------- */
484: /*
485:    Hessian Product - Computes the matrix-vector product:
486:    y = f''(x)*svec.

488:    Input Parameters
489: .  ptr  - pointer to the user-defined context
490: .  svec - input vector

492:    Output Parameters:
493: .  y    - product vector
494: */
495: PetscErrorCode HessianProduct(void *ptr,Vec svec,Vec y)
496: {
497:   AppCtx         *user = (AppCtx *)ptr;
498:   PetscReal      p5 = 0.5, zero = 0.0, one = 1.0, hx, hy, val, area, *x, *s;
499:   PetscReal      v, vb, vl, vr, vt, hxhx, hyhy;
501:   PetscInt       nx, ny, i, j, k, ind;

503:   nx   = user->mx;
504:   ny   = user->my;
505:   hx   = user->hx;
506:   hy   = user->hy;
507:   hxhx = one/(hx*hx);
508:   hyhy = one/(hy*hy);

510:   /* Get pointers to vector data */
511:   VecGetArray(user->xvec,&x);
512:   VecGetArray(svec,&s);

514:   /* Initialize product vector to zero */
515:   VecSet(y, zero);

517:   /* Compute f''(x)*s over the lower triangular elements */
518:   for (j=-1; j<ny; j++) {
519:     for (i=-1; i<nx; i++) {
520:        k = nx*j + i;
521:        v = zero;
522:        vr = zero;
523:        vt = zero;
524:        if (i != -1 && j != -1) v = s[k];
525:        if (i != nx-1 && j != -1) {
526:          vr = s[k+1];
527:          ind = k+1; val = hxhx*(vr-v);
528:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
529:        }
530:        if (i != -1 && j != ny-1) {
531:          vt = s[k+nx];
532:          ind = k+nx; val = hyhy*(vt-v);
533:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
534:        }
535:        if (i != -1 && j != -1) {
536:          ind = k; val = hxhx*(v-vr) + hyhy*(v-vt);
537:          VecSetValues(y,1,&ind,&val,ADD_VALUES);
538:        }
539:      }
540:    }

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

570:   /* Assemble vector */
571:   VecAssemblyBegin(y);
572:   VecAssemblyEnd(y);

574:   /* Scale resulting vector by area */
575:   area = p5*hx*hy;
576:   VecScale(y, area);
577:   PetscLogFlops(nx*ny*18);
578:   return 0;
579: }