Actual source code: jbearing2.c

      1: /*$Id: jbearing2.c 1.78 05/05/12 13:59:03-05:00 sarich@zorak.(none) $*/

      3: /*
      4:   Include "tao.h" so we can use TAO solvers with PETSc support.
      5:   Include "petscda.h" so that we can use distributed arrays (DAs) for managing
      6:   the parallel mesh.
      7: */

      9: #include "petscda.h"
     10:  #include tao.h
     11: #include <math.h>  /* for cos() sin(0), and atan() */

     13: static  char help[]=
     14: "This example demonstrates use of the TAO package to \n\
     15: solve a bound constrained minimization problem.  This example is based on \n\
     16: the problem DPJB from the MINPACK-2 test suite.  This pressure journal \n\
     17: bearing problem is an example of elliptic variational problem defined over \n\
     18: a two dimensional rectangle.  By discretizing the domain into triangular \n\
     19: elements, the pressure surrounding the journal bearing is defined as the \n\
     20: minimum of a quadratic function whose variables are bounded below by zero.\n\
     21: The command line options are:\n\
     22:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
     23:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
     24:  \n";

     26: /*T
     27:    Concepts: TAO - Solving a bound constrained minimization problem
     28:    Routines: TaoInitialize(); TaoFinalize();
     29:    Routines: TaoApplicationCreate();  TaoAppDestroy();
     30:    Routines: TaoAppSetObjectiveAndGradientRoutine(); TaoAppSetHessianRoutine();
     31:    Routines: TaoAppSetInitialSolutionVec(); TaoAppSetHessianMat();
     32:    Routines: TaoCreate(); TaoDestroy();
     33:    Routines: TaoSetOptions(); TaoGetKSP()
     34:    Routines: TaoSolveApplication(); TaoGetTerminationReason();
     35:    Processors: n
     36: T*/

     38: /*
     39:    User-defined application context - contains data needed by the
     40:    application-provided call-back routines, FormFunctionGradient(),
     41:    FormHessian().
     42: */
     43: typedef struct {
     44:   /* problem parameters */
     45:   PetscReal      ecc;          /* test problem parameter */
     46:   PetscReal      b;            /* A dimension of journal bearing */
     47:   int         nx,ny;        /* discretization in x, y directions */

     49:   /* Working space */
     50:   DA          da;           /* distributed array data structure */
     51:   Mat         A;            /* Quadratic Objective term */
     52:   Vec         B;            /* Linear Objective term */
     53: } AppCtx;

     55: /* User-defined routines */
     56: static PetscReal p(PetscReal xi, PetscReal ecc);
     57: static int FormFunctionGradient(TAO_APPLICATION, Vec, double *,Vec,void *);
     58: static int FormHessian(TAO_APPLICATION,Vec,Mat *, Mat *, MatStructure *, void *);
     59: static int ComputeB(AppCtx*);
     60: static int ComputeBounds(TAO_APPLICATION, Vec, Vec, void *);

     64: int main( int argc, char **argv )
     65: {
     66:   int        info;               /* used to check for functions returning nonzeros */
     67:   int        Nx, Ny;             /* number of processors in x- and y- directions */
     68:   int        m, N;               /* number of local and global elements in vectors */
     69:   Vec        x;                  /* variables vector */
     70:   PetscTruth   flg;              /* A return variable when checking for user options */
     71:   TAO_SOLVER tao;                /* TAO_SOLVER solver context */
     72:   TaoMethod  method = "tao_gpcg";/* default minimization method */
     73:   TAO_APPLICATION jbearingapp;   /* The PETSc application */
     74:   ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
     75:   int nloc;                      /* The number of local elements */
     76:   int *ltog;                     /* mapping of local elements to global elements */
     77:   TaoTerminateReason reason;
     78:   AppCtx     user;               /* user-defined work context */
     79:   PetscReal     zero=0.0;           /* lower bound on all variables */
     80:   KSP      ksp;

     82:
     83:   /* Initialize PETSC and TAO */
     84:   PetscInitialize( &argc, &argv,(char *)0,help );
     85:   TaoInitialize( &argc, &argv,(char *)0,help );

     87:   /* Set the default values for the problem parameters */
     88:   user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;

     90:   /* Check for any command line arguments that override defaults */
     91:   info = PetscOptionsGetInt(TAO_NULL,"-mx",&user.nx,&flg); CHKERRQ(info);
     92:   info = PetscOptionsGetInt(TAO_NULL,"-my",&user.ny,&flg); CHKERRQ(info);
     93:   info = PetscOptionsGetReal(TAO_NULL,"-ecc",&user.ecc,&flg); CHKERRQ(info);
     94:   info = PetscOptionsGetReal(TAO_NULL,"-b",&user.b,&flg); CHKERRQ(info);

     97:   PetscPrintf(MPI_COMM_WORLD,"\n---- Journal Bearing Problem -----\n");
     98:   PetscPrintf(MPI_COMM_WORLD,"mx: %d,  my: %d,  ecc: %4.3f \n\n",
     99:               user.nx,user.ny,user.ecc);

    101:   /* Calculate any derived values from parameters */
    102:   N = user.nx*user.ny;

    104:   /* Let Petsc determine the grid division */
    105:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

    107:   /*
    108:      A two dimensional distributed array will help define this problem,
    109:      which derives from an elliptic PDE on two dimensional domain.  From
    110:      the distributed array, Create the vectors.
    111:   */
    112:   info = DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.nx,
    113:                     user.ny,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da); CHKERRQ(info);

    115:   /*
    116:      Extract global and local vectors from DA; the vector user.B is
    117:      used solely as work space for the evaluation of the function,
    118:      gradient, and Hessian.  Duplicate for remaining vectors that are
    119:      the same types.
    120:   */
    121:   info = DACreateGlobalVector(user.da,&x); CHKERRQ(info); /* Solution */
    122:   info = VecDuplicate(x,&user.B); CHKERRQ(info); /* Linear objective */

    125:   /*  Create matrix user.A to store quadratic, Create a local ordering scheme. */
    126:   info = VecGetLocalSize(x,&m); CHKERRQ(info);
    127:   info = MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,5,TAO_NULL,3,TAO_NULL,&user.A); CHKERRQ(info);
    128:
    129:   info = DAGetGlobalIndices(user.da,&nloc,&ltog); CHKERRQ(info);
    130:   info = ISLocalToGlobalMappingCreate(MPI_COMM_SELF,nloc,ltog,&isltog);
    131:   CHKERRQ(info);
    132:   info = MatSetLocalToGlobalMapping(user.A,isltog); CHKERRQ(info);
    133:   info = ISLocalToGlobalMappingDestroy(isltog); CHKERRQ(info);

    136:   /* User defined function -- compute linear term of quadratic */
    137:   info = ComputeB(&user); CHKERRQ(info);

    140:   /* The TAO code begins here */

    142:   /*
    143:      Create the optimization solver, Petsc application
    144:      Suitable methods: "tao_gpcg","tao_bqpip","tao_tron","tao_blmvm"
    145:   */
    146:   info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
    147:   info = TaoApplicationCreate(MPI_COMM_WORLD,&jbearingapp); CHKERRQ(info);

    149:   /* Set the initial vector */
    150:   info = VecSet(x, zero); CHKERRQ(info);
    151:   info = TaoAppSetInitialSolutionVec(jbearingapp,x); CHKERRQ(info);

    153:   /* Set the user function, gradient, hessian evaluation routines and data structures */
    154:   info = TaoAppSetObjectiveAndGradientRoutine(jbearingapp,FormFunctionGradient,(void*) &user);
    155:   CHKERRQ(info);
    156:
    157:   info = TaoAppSetHessianMat(jbearingapp,user.A,user.A); CHKERRQ(info);
    158:   info = TaoAppSetHessianRoutine(jbearingapp,FormHessian,(void*)&user);
    159:   CHKERRQ(info);

    161:   /* Set a routine that defines the bounds */
    162:   info = TaoAppSetVariableBoundsRoutine(jbearingapp,ComputeBounds,(void*)&user); CHKERRQ(info);

    164:   info = TaoGetKSP(tao,&ksp); CHKERRQ(info);
    165:   if (ksp) {                                              /* Modify the PETSc KSP structure */
    166:     info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
    167:   }

    170:   /* Check for any tao command line options */
    171:   info = TaoSetOptions(jbearingapp,tao); CHKERRQ(info);

    173:   /* Solve the bound constrained problem */
    174:   info = TaoSolveApplication(jbearingapp,tao); CHKERRQ(info);

    176:   info = TaoGetTerminationReason(tao,&reason); CHKERRQ(info);
    177:   if (reason <= 0)
    178:     PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");

    180:   /* Free TAO data structures */
    181:   info = TaoDestroy(tao); CHKERRQ(info);
    182:   info = TaoAppDestroy(jbearingapp); CHKERRQ(info);

    184:   /* Free PETSc data structures */
    185:   info = VecDestroy(x); CHKERRQ(info);
    186:   info = MatDestroy(user.A); CHKERRQ(info);
    187:   info = VecDestroy(user.B); CHKERRQ(info);
    188:   info = DADestroy(user.da); CHKERRQ(info);

    190:   TaoFinalize();
    191:   PetscFinalize();

    193:   return 0;
    194: }

    198: static int ComputeBounds(TAO_APPLICATION taoapp, Vec xl, Vec xu, void *ctx){
    199:   int info;
    200:   PetscReal zero=0.0, d1000=1000;
    201:   /* Set the upper and lower bounds */
    202:   info = VecSet(xl, zero); CHKERRQ(info);
    203:   info = VecSet(xu, d1000); CHKERRQ(info);
    204:   return 0;
    205: }

    207: static PetscReal p(PetscReal xi, PetscReal ecc)
    208: {
    209:   PetscReal t=1.0+ecc*cos(xi);
    210:   return (t*t*t);
    211: }

    215: int ComputeB(AppCtx* user)
    216: {
    217:   int i,j,k,info;
    218:   int nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
    219:   PetscReal two=2.0, pi=4.0*atan(1.0);
    220:   PetscReal hx,hy,ehxhy;
    221:   PetscReal temp,*b;
    222:   PetscReal ecc=user->ecc;

    224:   nx=user->nx;
    225:   ny=user->ny;
    226:   hx=two*pi/(nx+1.0);
    227:   hy=two*user->b/(ny+1.0);
    228:   ehxhy = ecc*hx*hy;

    231:   /*
    232:      Get local grid boundaries
    233:   */
    234:   info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
    235:   info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
    236:

    238:   /* Compute the linear term in the objective function */
    239:   info = VecGetArray(user->B,&b); CHKERRQ(info);
    240:   for (i=xs; i<xs+xm; i++){
    241:     temp=sin((i+1)*hx);
    242:     for (j=ys; j<ys+ym; j++){
    243:       k=xm*(j-ys)+(i-xs);
    244:       b[k]=  - ehxhy*temp;
    245:     }
    246:   }
    247:   info = VecRestoreArray(user->B,&b); CHKERRQ(info);
    248:   info = PetscLogFlops(5*xm*ym+3*xm); CHKERRQ(info);

    250:   return 0;
    251: }

    255: int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn,Vec G,void *ptr)
    256: {
    257:   AppCtx* user=(AppCtx*)ptr;
    258:   int i,j,k,kk,info;
    259:   int col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
    260:   PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
    261:   PetscReal hx,hy,hxhy,hxhx,hyhy;
    262:   PetscReal xi,v[5];
    263:   PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
    264:   PetscReal vmiddle, vup, vdown, vleft, vright;
    265:   PetscReal tt,f1,f2;
    266:   PetscReal *x,*g,zero=0.0;
    267:   Vec localX;

    269:   nx=user->nx;
    270:   ny=user->ny;
    271:   hx=two*pi/(nx+1.0);
    272:   hy=two*user->b/(ny+1.0);
    273:   hxhy=hx*hy;
    274:   hxhx=one/(hx*hx);
    275:   hyhy=one/(hy*hy);

    277:   info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);

    279:   info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
    280:   info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);

    282:   info = VecSet(G, zero); CHKERRQ(info);
    283:   /*
    284:     Get local grid boundaries
    285:   */
    286:   info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
    287:   info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
    288:
    289:   info = VecGetArray(localX,&x); CHKERRQ(info);
    290:   info = VecGetArray(G,&g); CHKERRQ(info);

    292:   for (i=xs; i< xs+xm; i++){
    293:     xi=(i+1)*hx;
    294:     trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
    295:     trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
    296:     trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
    297:     trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
    298:     trule5=trule1; /* L(i,j-1) */
    299:     trule6=trule2; /* U(i,j+1) */

    301:     vdown=-(trule5+trule2)*hyhy;
    302:     vleft=-hxhx*(trule2+trule4);
    303:     vright= -hxhx*(trule1+trule3);
    304:     vup=-hyhy*(trule1+trule6);
    305:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);

    307:     for (j=ys; j<ys+ym; j++){
    308:
    309:       row=(j-gys)*gxm + (i-gxs);
    310:        v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
    311:
    312:        k=0;
    313:        if (j>gys){
    314:          v[k]=vdown; col[k]=row - gxm; k++;
    315:        }
    316:
    317:        if (i>gxs){
    318:          v[k]= vleft; col[k]=row - 1; k++;
    319:        }

    321:        v[k]= vmiddle; col[k]=row; k++;
    322:
    323:        if (i+1 < gxs+gxm){
    324:          v[k]= vright; col[k]=row+1; k++;
    325:        }
    326:
    327:        if (j+1 <gys+gym){
    328:          v[k]= vup; col[k] = row+gxm; k++;
    329:        }
    330:        tt=0;
    331:        for (kk=0;kk<k;kk++){
    332:          tt+=v[kk]*x[col[kk]];
    333:        }
    334:        row=(j-ys)*xm + (i-xs);
    335:        g[row]=tt;

    337:      }

    339:   }

    341:   info = VecRestoreArray(localX,&x); CHKERRQ(info);
    342:   info = VecRestoreArray(G,&g); CHKERRQ(info);

    344:   info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);

    346:   info = VecDot(X,G,&f1); CHKERRQ(info);
    347:   info = VecDot(user->B,X,&f2); CHKERRQ(info);
    348:   info = VecAXPY(G, one, user->B); CHKERRQ(info);
    349:   *fcn = f1/2.0 + f2;

    351:   info = PetscLogFlops((91 + 10*ym) * xm); CHKERRQ(info);
    352:   return 0;

    354: }

    360: /*
    361:    FormHessian computes the quadratic term in the quadratic objective function
    362:    Notice that the objective function in this problem is quadratic (therefore a constant
    363:    hessian).  If using a nonquadratic solver, then you might want to reconsider this function
    364: */
    365: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
    366: {
    367:   AppCtx* user=(AppCtx*)ptr;
    368:   int i,j,k,info;
    369:   int col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
    370:   PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
    371:   PetscReal hx,hy,hxhy,hxhx,hyhy;
    372:   PetscReal xi,v[5];
    373:   PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
    374:   PetscReal vmiddle, vup, vdown, vleft, vright;
    375:   Mat hes=*H;
    376:   PetscTruth assembled;

    378:   nx=user->nx;
    379:   ny=user->ny;
    380:   hx=two*pi/(nx+1.0);
    381:   hy=two*user->b/(ny+1.0);
    382:   hxhy=hx*hy;
    383:   hxhx=one/(hx*hx);
    384:   hyhy=one/(hy*hy);

    386:   *flg=SAME_NONZERO_PATTERN;
    387:   /*
    388:     Get local grid boundaries
    389:   */
    390:   info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
    391:   info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
    392:
    393:   info = MatAssembled(hes,&assembled); CHKERRQ(info);
    394:   if (assembled){info = MatZeroEntries(hes);  CHKERRQ(info);}

    396:   for (i=xs; i< xs+xm; i++){
    397:     xi=(i+1)*hx;
    398:     trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
    399:     trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
    400:     trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
    401:     trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
    402:     trule5=trule1; /* L(i,j-1) */
    403:     trule6=trule2; /* U(i,j+1) */

    405:     vdown=-(trule5+trule2)*hyhy;
    406:     vleft=-hxhx*(trule2+trule4);
    407:     vright= -hxhx*(trule1+trule3);
    408:     vup=-hyhy*(trule1+trule6);
    409:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
    410:     v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;

    412:     for (j=ys; j<ys+ym; j++){
    413:       row=(j-gys)*gxm + (i-gxs);
    414:
    415:       k=0;
    416:       if (j>gys){
    417:         v[k]=vdown; col[k]=row - gxm; k++;
    418:       }
    419:
    420:       if (i>gxs){
    421:         v[k]= vleft; col[k]=row - 1; k++;
    422:       }

    424:       v[k]= vmiddle; col[k]=row; k++;
    425:
    426:       if (i+1 < gxs+gxm){
    427:         v[k]= vright; col[k]=row+1; k++;
    428:       }
    429:
    430:       if (j+1 <gys+gym){
    431:         v[k]= vup; col[k] = row+gxm; k++;
    432:       }
    433:       info = MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info);
    434:
    435:     }

    437:   }

    439:   /*
    440:      Assemble matrix, using the 2-step process:
    441:      MatAssemblyBegin(), MatAssemblyEnd().
    442:      By placing code between these two statements, computations can be
    443:      done while messages are in transition.
    444:   */
    445:   info = MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
    446:   info = MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY); CHKERRQ(info);

    448:   /*
    449:     Tell the matrix we will never add a new nonzero location to the
    450:     matrix. If we do it will generate an error.
    451:   */
    452:   info = MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ(info);
    453:   info = MatSetOption(hes,MAT_SYMMETRIC); CHKERRQ(info);

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

    457:   return 0;
    458: }