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();
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
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 */
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: /*
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: }
```