Actual source code: eptorsion2.c

      1: /*$Id: eptorsion2.c 1.55 05/05/11 08:44:55-05:00 sarich@zorak.(none) $*/

      3: /* Program usage: mpirun -np <proc> eptorsion2 [-help] [all TAO options] */

      5: /* ----------------------------------------------------------------------

      7:   Elastic-plastic torsion problem.

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

     13:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
     14:
     15:   where C is the torsion angle per unit length.

     17:   The uniprocessor version of this code is eptorsion1.c; the Fortran
     18:   version of this code is eptorsion2f.F.

     20:   This application solves the problem without calculating hessians
     21: ---------------------------------------------------------------------- */

     23: /*
     24:   Include "tao.h" so that we can use TAO solvers.  Note that this
     25:   file automatically includes files for lower-level support, such as those
     26:   provided by the PETSc library:
     27:      petsc.h       - base PETSc routines   petscvec.h - vectors
     28:      petscsys.h    - sysem routines        petscmat.h - matrices
     29:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     30:      petscviewer.h - viewers               petscpc.h  - preconditioners
     31:   Include "petscda.h" so that we can use distributed arrays (DAs) for managing
     32:   the parallel mesh.
     33: */

     35:  #include tao.h
     36: #include "petscda.h"

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

     47: /*T
     48:    Concepts: TAO - Solving an unconstrained minimization problem
     49:    Routines: TaoInitialize(); TaoFinalize();
     50:    Routines: TaoApplicationCreate(); TaoAppDestroy();
     51:    Routines: TaoCreate(); TaoDestroy();
     52:    Routines: TaoAppSetObjectiveAndGradientRoutine();
     53:    Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
     54:    Routines: TaoSetOptions();
     55:    Routines: TaoAppSetInitialSolutionVec();
     56:    Routines: TaoSolve(); TaoView(); TaoGetKSP();
     57:    Routines: TaoGetSolutionStatus();
     58:    Processors: n
     59: T*/

     63: /*
     64:    User-defined application context - contains data needed by the
     65:    application-provided call-back routines, FormFunction() and
     66:    FormGradient().
     67: */
     68: typedef struct {
     69:   /* parameters */
     70:    int           mx, my;         /* global discretization in x- and y-directions */
     71:    PetscReal        param;          /* nonlinearity parameter */

     73:   /* work space */
     74:    Vec           localX;         /* local vectors */
     75:    DA            da;             /* distributed array data structure */
     76: } AppCtx;

     78: /* -------- User-defined Routines --------- */

     80: int FormInitialGuess(AppCtx*,Vec);
     81: int FormFunctionGradient(TAO_APPLICATION,Vec,double*,Vec,void*);
     82: int ComputeHessian(TAO_APPLICATION,Vec,Mat*,Mat*,MatStructure*,void*);

     86: int main(int argc,char **argv)
     87: {
     88:   int             info;                  /* used to check for functions returning nonzeros */
     89:   Vec             x;                     /* solution, gradient vectors */
     90:   Mat             H;                     /* Hessian matrix */
     91:   int             Nx, Ny;                /* number of processes in x- and y-directions */
     92:   TAO_SOLVER      tao;                   /* TAO_SOLVER solver context */
     93:   TAO_APPLICATION torsionapp;            /* TAO application context */
     94:   TaoTerminateReason reason;
     95:   PetscTruth      flg;
     96:   int             iter;                  /* iteration information */
     97:   double          ff,gnorm;
     98:   AppCtx          user;                  /* application context */
     99:   KSP             ksp;

    101:   /* Initialize TAO, PETSc contexts */
    102:   info = PetscInitialize(&argc,&argv,(char *)0,help);
    103:   info = TaoInitialize(&argc,&argv,(char *)0,help);

    105:   /* Specify default parameters */
    106:   user.param = 5.0; user.mx = 10; user.my = 10;
    107:   Nx = Ny = PETSC_DECIDE;

    109:   /* Check for any command line arguments that override defaults */
    110:   info = PetscOptionsGetReal(TAO_NULL,"-par",&user.param,&flg); CHKERRQ(info);
    111:   info = PetscOptionsGetInt(TAO_NULL,"-my",&user.my,&flg); CHKERRQ(info);
    112:   info = PetscOptionsGetInt(TAO_NULL,"-mx",&user.mx,&flg); CHKERRQ(info);

    114:   /* Set up distributed array */
    115:   info = DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,
    116:                     user.mx,user.my,Nx,Ny,1,1,TAO_NULL,TAO_NULL,
    117:                     &user.da); CHKERRQ(info);

    119:   /* Create vectors */
    120:   info = DACreateGlobalVector(user.da,&x); CHKERRQ(info);

    122:   info = DACreateLocalVector(user.da,&user.localX); CHKERRQ(info);

    124:   /* Create Hessian */
    125:   info = DAGetMatrix(user.da,MATAIJ,&H); CHKERRQ(info);
    126:   info = MatSetOption(H,MAT_SYMMETRIC); CHKERRQ(info);

    128:   /* The TAO code begins here */

    130:   /* Create TAO solver and set desired solution method */
    131:   info = TaoCreate(MPI_COMM_WORLD,"tao_cg_fr",&tao); CHKERRQ(info);
    132:   info = TaoApplicationCreate(PETSC_COMM_WORLD,&torsionapp); CHKERRQ(info);

    134:   /* Set initial solution guess */
    135:   info = FormInitialGuess(&user,x); CHKERRQ(info);
    136:   info = TaoAppSetInitialSolutionVec(torsionapp,x); CHKERRQ(info);

    138:   /* Set routine for function and gradient evaluation */
    139:   info = TaoAppSetObjectiveAndGradientRoutine(torsionapp,FormFunctionGradient,(void *)&user); CHKERRQ(info);

    141:   info = TaoAppSetHessianMat(torsionapp,H,H); CHKERRQ(info);
    142:   info = TaoAppSetHessianRoutine(torsionapp,ComputeHessian,(void*)&user); CHKERRQ(info);

    144:   info = TaoGetKSP(tao,&ksp); CHKERRQ(info);
    145:   if (ksp) {                                              /* Modify the PETSc KSP structure */
    146:     info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
    147:   }

    149:   /* Check for any TAO command line options */
    150:   info = TaoSetOptions(torsionapp,tao); CHKERRQ(info);

    152:   /* SOLVE THE APPLICATION */
    153:   info = TaoSolveApplication(torsionapp,tao);  CHKERRQ(info);

    155:   /* Get information on termination */
    156:   info = TaoGetSolutionStatus(tao,&iter,&ff,&gnorm,0,0,&reason); CHKERRQ(info);
    157:   if (reason <= 0){
    158:     PetscPrintf(MPI_COMM_WORLD, "Try another method! Iterations: %d, f: %4.2e, residual: %4.2e\n",
    159:                 iter,ff,gnorm); CHKERRQ(info);
    160:   }

    162:   /*
    163:      To View TAO solver information use
    164:       info = TaoView(tao); CHKERRQ(info);
    165:   */

    167:   /* Free TAO data structures */
    168:   info = TaoDestroy(tao); CHKERRQ(info);
    169:   info = TaoAppDestroy(torsionapp);  CHKERRQ(info);

    171:   /* Free PETSc data structures */
    172:   info = VecDestroy(x); CHKERRQ(info);
    173:   info = MatDestroy(H); CHKERRQ(info);

    175:   info = VecDestroy(user.localX); CHKERRQ(info);
    176:   info = DADestroy(user.da); CHKERRQ(info);

    179:   /* Finalize TAO and PETSc*/
    180:   TaoFinalize();
    181:   PetscFinalize();

    183:   return 0;
    184: }

    187: /* ------------------------------------------------------------------- */
    190: /*
    191:     FormInitialGuess - Computes an initial approximation to the solution.

    193:     Input Parameters:
    194: .   user - user-defined application context
    195: .   X    - vector

    197:     Output Parameters:
    198:     X    - vector
    199: */
    200: int FormInitialGuess(AppCtx *user,Vec X)
    201: {
    202:   int    info, i, j, k, mx = user->mx, my = user->my;
    203:   int    xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
    204:   PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;

    206:   /* Get local mesh boundaries */
    207:   info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
    208:   info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);

    210:   /* Compute initial guess over locally owned part of mesh */
    211:   xe = xs+xm;
    212:   ye = ys+ym;
    213:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
    214:     temp = TaoMin(j+1,my-j)*hy;
    215:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
    216:       k   = (j-gys)*gxm + i-gxs;
    217:       val = TaoMin((TaoMin(i+1,mx-i))*hx,temp);
    218:       info = VecSetValuesLocal(X,1,&k,&val,ADD_VALUES); CHKERRQ(info);
    219:     }
    220:   }

    222:   return 0;
    223: }

    226: /* ------------------------------------------------------------------ */
    229: /*
    230:    FormFunctionGradient - Evaluates the function and corresponding gradient.
    231:
    232:    Input Parameters:
    233:    taoapp - the TAO_APPLICATION context
    234:    X   - the input vector
    235:    ptr - optional user-defined context, as set by TaoSetFunction()

    237:    Output Parameters:
    238:    f   - the newly evaluated function
    239:    G   - the newly evaluated gradient
    240: */
    241: int FormFunctionGradient(TAO_APPLICATION taoapp,Vec X,double *f,Vec G,void *ptr){

    243:   AppCtx *user = (AppCtx *)ptr;
    244:   int    info,i,j,k,ind;
    245:   int    xe,ye,xsm,ysm,xep,yep;
    246:   int    xs, ys, xm, ym, gxm, gym, gxs, gys;
    247:   int    mx = user->mx, my = user->my;
    248:   PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
    249:   PetscReal p5 = 0.5, area, val, flin, fquad;
    250:   PetscReal v,vb,vl,vr,vt,dvdx,dvdy;
    251:   PetscReal hx = 1.0/(user->mx + 1);
    252:   PetscReal hy = 1.0/(user->my + 1);
    253:   Vec    localX = user->localX;

    256:   /* Initialize */
    257:   flin = fquad = zero;

    259:   info = VecSet(G, zero); CHKERRQ(info);
    260:   /*
    261:      Scatter ghost points to local vector,using the 2-step process
    262:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
    263:      By placing code between these two statements, computations can be
    264:      done while messages are in transition.
    265:   */
    266:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
    267:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);

    269:   /* Get pointer to vector data */
    270:   info = VecGetArray(localX,&x); CHKERRQ(info);

    272:   /* Get local mesh boundaries */
    273:   info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
    274:   info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);

    276:   /* Set local loop dimensions */
    277:   xe = xs+xm;
    278:   ye = ys+ym;
    279:   if (xs == 0)  xsm = xs-1;
    280:   else          xsm = xs;
    281:   if (ys == 0)  ysm = ys-1;
    282:   else          ysm = ys;
    283:   if (xe == mx) xep = xe+1;
    284:   else          xep = xe;
    285:   if (ye == my) yep = ye+1;
    286:   else          yep = ye;

    288:   /* Compute local gradient contributions over the lower triangular elements */
    289:   for (j=ysm; j<ye; j++) {  /*  for (j=-1; j<my; j++) */
    290:     for (i=xsm; i<xe; i++) {  /*   for (i=-1; i<mx; i++) */
    291:       k = (j-gys)*gxm + i-gxs;
    292:       v = zero;
    293:       vr = zero;
    294:       vt = zero;
    295:       if (i >= 0 && j >= 0) v = x[k];
    296:       if (i < mx-1 && j > -1) vr = x[k+1];
    297:       if (i > -1 && j < my-1) vt = x[k+gxm];
    298:       dvdx = (vr-v)/hx;
    299:       dvdy = (vt-v)/hy;
    300:       if (i != -1 && j != -1) {
    301:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
    302:         info = VecSetValuesLocal(G,1,&k,&val,ADD_VALUES); CHKERRQ(info);
    303:       }
    304:       if (i != mx-1 && j != -1) {
    305:         ind = k+1; val =  dvdx/hx - cdiv3;
    306:         info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    307:       }
    308:       if (i != -1 && j != my-1) {
    309:         ind = k+gxm; val = dvdy/hy - cdiv3;
    310:         info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    311:       }
    312:       fquad += dvdx*dvdx + dvdy*dvdy;
    313:       flin -= cdiv3 * (v + vr + vt);
    314:     }
    315:   }

    317:   /* Compute local gradient contributions over the upper triangular elements */
    318:   for (j=ys; j<yep; j++) { /*  for (j=0; j<=my; j++) */
    319:     for (i=xs; i<xep; i++) {  /*   for (i=0; i<=mx; i++) */
    320:       k = (j-gys)*gxm + i-gxs;
    321:       vb = zero;
    322:       vl = zero;
    323:       v  = zero;
    324:       if (i < mx && j > 0) vb = x[k-gxm];
    325:       if (i > 0 && j < my) vl = x[k-1];
    326:       if (i < mx && j < my) v = x[k];
    327:       dvdx = (v-vl)/hx;
    328:       dvdy = (v-vb)/hy;
    329:       if (i != mx && j != 0) {
    330:         ind = k-gxm; val = - dvdy/hy - cdiv3;
    331:         info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    332:       }
    333:       if (i != 0 && j != my) {
    334:         ind = k-1; val =  - dvdx/hx - cdiv3;
    335:         info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    336:       }
    337:       if (i != mx && j != my) {
    338:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
    339:         info = VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); CHKERRQ(info);
    340:       }
    341:       fquad += dvdx*dvdx + dvdy*dvdy;
    342:       flin -= cdiv3 * (vb + vl + v);
    343:     }
    344:   }

    347:   /* Restore vector */
    348:   info = VecRestoreArray(localX,&x); CHKERRQ(info);

    350:   /* Assemble gradient vector */
    351:   info = VecAssemblyBegin(G); CHKERRQ(info);
    352:   info = VecAssemblyEnd(G); CHKERRQ(info);

    354:   /* Scale the gradient */
    355:   area = p5*hx*hy;
    356:   floc = area * (p5 * fquad + flin);
    357:   info = VecScale(G, area); CHKERRQ(info);

    359:   /* Sum function contributions from all processes */
    360:   MPI_Allreduce((void*)&floc,(void*)f,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);

    362:   info=PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16); CHKERRQ(info);
    363:
    364:   return 0;
    365: }

    371: int ComputeHessian(TAO_APPLICATION taoapp, Vec X, Mat *H, Mat *Hpre, MatStructure *flag, void*ctx){

    373:   AppCtx *user= (AppCtx*) ctx;
    374:   int i,j,k,info;
    375:   int col[5],row;
    376:   int xs,xm,gxs,gxm,ys,ym,gys,gym;
    377:   PetscReal v[5];
    378:   PetscReal hx=1.0/(user->mx+1), hy=1.0/(user->my+1), hxhx=1.0/(hx*hx), hyhy=1.0/(hy*hy), area=0.5*hx*hy;
    379:   Mat A=*H;

    381:   /* Compute the quadratic term in the objective function */

    383:   /*
    384:      Get local grid boundaries
    385:   */

    387:   info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
    388:   info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);

    390:   for (j=ys; j<ys+ym; j++){
    391:
    392:     for (i=xs; i< xs+xm; i++){

    394:       row=(j-gys)*gxm + (i-gxs);

    396:       k=0;
    397:       if (j>gys){
    398:         v[k]=-2*hyhy; col[k]=row - gxm; k++;
    399:       }

    401:       if (i>gxs){
    402:         v[k]= -2*hxhx; col[k]=row - 1; k++;
    403:       }

    405:       v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++;

    407:       if (i+1 < gxs+gxm){
    408:         v[k]= -2.0*hxhx; col[k]=row+1; k++;
    409:       }

    411:       if (j+1 <gys+gym){
    412:         v[k]= -2*hyhy; col[k] = row+gxm; k++;
    413:       }

    415:       info = MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info);

    417:     }
    418:   }
    419:   /*
    420:      Assemble matrix, using the 2-step process:
    421:      MatAssemblyBegin(), MatAssemblyEnd().
    422:      By placing code between these two statements, computations can be
    423:      done while messages are in transition.
    424:   */
    425:   info = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    426:   info = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    427:   /*
    428:     Tell the matrix we will never add a new nonzero location to the
    429:     matrix. If we do it will generate an error.
    430:   */
    431:   info = MatScale(A,area); CHKERRQ(info);
    432:   info = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ(info);
    433:   info = MatSetOption(A,MAT_SYMMETRIC); CHKERRQ(info);

    435:   info = PetscLogFlops(9*xm*ym+49*xm); CHKERRQ(info);

    437:   return 0;
    438: }