Actual source code: jbearing2.c

tao-2.1-p0 2012-07-24
  1: /*
  2:   Include "taosolver.h" so we can use TAO solvers
  3:   Include "petscdm.h" so that we can use distributed arrays (DMs) for managing
  4:   Include "petscksp.h" so we can set KSP type
  5:   the parallel mesh.
  6: */

  8: #include "petscdm.h"
  9: #include "petscksp.h"
 10: #include "taosolver.h"

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

 25: /*T
 26:    Concepts: TAO - Solving a bound constrained minimization problem
 27:    Routines: TaoInitialize(); TaoFinalize();
 28:    Routines: TaoCreate();
 29:    Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine();
 30:    Routines: TaoSetHessianRoutine();
 31:    Routines: TaoSetVariableBounds();
 32:    Routines: TaoSetMonitor(); TaoSetConvergenceTest();
 33:    Routines: TaoSetInitialVector();
 34:    Routines: TaoSetFromOptions();
 35:    Routines: TaoSolve();
 36:    Routines: TaoGetTerminationReason(); TaoDestroy(); 
 37:    Processors: n
 38: T*/

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

 51:   /* Working space */
 52:   DM          dm;           /* distributed array data structure */
 53:   Mat         A;            /* Quadratic Objective term */
 54:   Vec         B;            /* Linear Objective term */
 55: } AppCtx;

 57: /* User-defined routines */
 58: static PetscReal p(PetscReal xi, PetscReal ecc);
 59: static PetscErrorCode FormFunctionGradient(TaoSolver, Vec, PetscReal *,Vec,void *);
 60: static PetscErrorCode FormHessian(TaoSolver,Vec,Mat *, Mat *, MatStructure *, void *);
 61: static PetscErrorCode ComputeB(AppCtx*);
 62: static PetscErrorCode Monitor(TaoSolver, void*);
 63: static PetscErrorCode ConvergenceTest(TaoSolver, void*);

 67: int main( int argc, char **argv )
 68: {
 69:   PetscErrorCode        ierr;               /* used to check for functions returning nonzeros */
 70:   PetscInt        Nx, Ny;             /* number of processors in x- and y- directions */
 71:   PetscInt        m;               /* number of local elements in vectors */
 72:   Vec        x;                  /* variables vector */
 73:   Vec        xl,xu;                  /* bounds vectors */
 74:   PetscReal d1000 = 1000;
 75:   PetscBool   flg;              /* A return variable when checking for user options */
 76:   TaoSolver tao;                /* TaoSolver solver context */

 78:   TaoSolverTerminationReason reason;
 79:   KSP        ksp;
 80:   AppCtx     user;               /* user-defined work context */
 81:   PetscReal     zero=0.0;           /* lower bound on all variables */


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

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

 92:   /* Check for any command line arguments that override defaults */
 93:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.nx,&flg); 
 94:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.ny,&flg); 
 95:   PetscOptionsGetReal(PETSC_NULL,"-ecc",&user.ecc,&flg); 
 96:   PetscOptionsGetReal(PETSC_NULL,"-b",&user.b,&flg); 


 99:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Journal Bearing Problem SHB-----\n");
100:   PetscPrintf(PETSC_COMM_WORLD,"mx: %D,  my: %D,  ecc: %G \n\n",
101:               user.nx,user.ny,user.ecc);

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

106:   /*
107:      A two dimensional distributed array will help define this problem,
108:      which derives from an elliptic PDE on two dimensional domain.  From
109:      the distributed array, Create the vectors.
110:   */
111:   DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
112:                       user.nx,user.ny,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,
113:                       &user.dm); 

115:   /*
116:      Extract global and local vectors from DM; 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:   DMCreateGlobalVector(user.dm,&x);  /* Solution */
122:   VecDuplicate(x,&user.B);  /* Linear objective */


125:   /*  Create matrix user.A to store quadratic, Create a local ordering scheme. */
126:   VecGetLocalSize(x,&m); 
127:   DMCreateMatrix(user.dm,MATAIJ,&user.A); 

129:   /* User defined function -- compute linear term of quadratic */
130:   ComputeB(&user); 

132:   /* The TAO code begins here */

134:   /* 
135:      Create the optimization solver, Petsc application 
136:      Suitable methods: "tao_gpcg","tao_bqpip","tao_tron","tao_blmvm" 
137:   */
138:   TaoCreate(PETSC_COMM_WORLD,&tao); 
139:   TaoSetType(tao,"tao_blmvm"); 


142:   /* Set the initial vector */
143:   VecSet(x, zero); 
144:   TaoSetInitialVector(tao,x); 

146:   /* Set the user function, gradient, hessian evaluation routines and data structures */
147:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user); 
148:   
149:   
150:   TaoSetHessianRoutine(tao,user.A,user.A,FormHessian,(void*)&user); 

152:   /* Set a routine that defines the bounds */
153:   VecDuplicate(x,&xl); 
154:   VecDuplicate(x,&xu); 
155:   VecSet(xl, zero); 
156:   VecSet(xu, d1000); 
157:   TaoSetVariableBounds(tao,xl,xu); 

159:   TaoGetKSP(tao,&ksp); 
160:   if (ksp) {                                         
161:     KSPSetType(ksp,KSPCG); 
162:   }

164:   PetscOptionsHasName(PETSC_NULL,"-testmonitor",&flg); 
165:   if (flg) {
166:     TaoSetMonitor(tao,Monitor,&user,PETSC_NULL); 
167:   }
168:   PetscOptionsHasName(PETSC_NULL,"-testconvergence",&flg);
169:   if (flg) {
170:     TaoSetConvergenceTest(tao,ConvergenceTest,&user); 
171:   }

173:   /* Check for any tao command line options */
174:   TaoSetFromOptions(tao); 

176:   /* Solve the bound constrained problem */
177:   TaoSolve(tao); 

179:   TaoGetTerminationReason(tao,&reason); 
180:   if (reason <= 0)
181:     PetscPrintf(PETSC_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");


184:   /* Free PETSc data structures */
185:   VecDestroy(&x);  
186:   VecDestroy(&xl);  
187:   VecDestroy(&xu);  
188:   MatDestroy(&user.A); 
189:   VecDestroy(&user.B);  
190:   /* Free TAO data structures */
191:   TaoDestroy(&tao); 

193:   DMDestroy(&user.dm); 

195:   TaoFinalize();
196:   PetscFinalize();

198:   return 0;
199: }


202: static PetscReal p(PetscReal xi, PetscReal ecc)
203: { 
204:   PetscReal t=1.0+ecc*PetscCosScalar(xi); 
205:   return (t*t*t); 
206: }

210: PetscErrorCode ComputeB(AppCtx* user)
211: {
213:   PetscInt i,j,k;
214:   PetscInt nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
215:   PetscReal two=2.0, pi=4.0*atan(1.0);
216:   PetscReal hx,hy,ehxhy;
217:   PetscReal temp,*b;
218:   PetscReal ecc=user->ecc;

220:   nx=user->nx;
221:   ny=user->ny;
222:   hx=two*pi/(nx+1.0);
223:   hy=two*user->b/(ny+1.0);
224:   ehxhy = ecc*hx*hy;


227:   /*
228:      Get local grid boundaries
229:   */
230:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
231:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 
232:   

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

246:   return 0;
247: }

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

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

274:   DMGetLocalVector(user->dm,&localX);

276:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
277:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

279:   VecSet(G, zero); 
280:   /*
281:     Get local grid boundaries
282:   */
283:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
284:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 
285:   
286:   VecGetArray(localX,&x); 
287:   VecGetArray(G,&g); 

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

298:     vdown=-(trule5+trule2)*hyhy;
299:     vleft=-hxhx*(trule2+trule4);
300:     vright= -hxhx*(trule1+trule3);
301:     vup=-hyhy*(trule1+trule6);
302:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);

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

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

334:      }

336:   }

338:   VecRestoreArray(localX,&x); 
339:   VecRestoreArray(G,&g); 

341:   DMRestoreLocalVector(user->dm,&localX); 

343:   VecDot(X,G,&f1); 
344:   VecDot(user->B,X,&f2); 
345:   VecAXPY(G, one, user->B); 
346:   *fcn = f1/2.0 + f2;
347:   

349:   PetscLogFlops((91 + 10*ym) * xm); 
350:   return 0;

352: }


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

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

384:   *flg=SAME_NONZERO_PATTERN;
385:   /*
386:     Get local grid boundaries
387:   */
388:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
389:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 
390:   
391:   MatAssembled(hes,&assembled); 
392:   if (assembled){MatZeroEntries(hes);  }

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

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

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

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

435:   }

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

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

453:   PetscLogFlops(9*xm*ym+49*xm); 
454:   MatNorm(hes,NORM_1,&hx); 
455:   return 0;
456: }


461: PetscErrorCode Monitor(TaoSolver tao, void *ctx)
462: {
464:   PetscInt its;
465:   PetscReal f,gnorm,cnorm,xdiff;
466:   TaoSolverTerminationReason reason;
468:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason); 
469:   if (!(its%5)) {
470:     PetscPrintf(PETSC_COMM_WORLD,"iteration=%D\tf=%G\n",its,f);
471:   }
472:   return(0);
473: }

477: PetscErrorCode ConvergenceTest(TaoSolver tao, void *ctx)
478: {
480:   PetscInt its;
481:   PetscReal f,gnorm,cnorm,xdiff;
482:   TaoSolverTerminationReason reason;
484:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason); 
485:   if (its == 100) {
486:     TaoSetTerminationReason(tao,TAO_DIVERGED_MAXITS);
487:   }
488:   return(0);
489:   
490: }