Actual source code: jbearing2.c

  1: /*
  2:   Include "petsctao.h" so we can use TAO solvers
  3:   Include "petscdmda.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 <petsctao.h>
  9: #include <petscdmda.h>

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

 24: /*T
 25:    Concepts: TAO^Solving a bound constrained minimization problem
 26:    Routines: TaoCreate();
 27:    Routines: TaoSetType(); TaoSetObjectiveAndGradient();
 28:    Routines: TaoSetHessian();
 29:    Routines: TaoSetVariableBounds();
 30:    Routines: TaoSetMonitor(); TaoSetConvergenceTest();
 31:    Routines: TaoSetSolution();
 32:    Routines: TaoSetFromOptions();
 33:    Routines: TaoSolve();
 34:    Routines: TaoDestroy();
 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:   PetscInt       nx,ny;        /* discretization in x, y directions */

 49:   /* Working space */
 50:   DM          dm;           /* 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 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *,Vec,void *);
 58: static PetscErrorCode FormHessian(Tao,Vec,Mat, Mat, void *);
 59: static PetscErrorCode ComputeB(AppCtx*);
 60: static PetscErrorCode Monitor(Tao, void*);
 61: static PetscErrorCode ConvergenceTest(Tao, void*);

 63: int main(int argc, char **argv)
 64: {
 65:   PetscInt           Nx, Ny;          /* number of processors in x- and y- directions */
 66:   PetscInt           m;               /* number of local elements in vectors */
 67:   Vec                x;               /* variables vector */
 68:   Vec                xl,xu;           /* bounds vectors */
 69:   PetscReal          d1000 = 1000;
 70:   PetscBool          flg,testgetdiag; /* A return variable when checking for user options */
 71:   Tao                tao;             /* Tao solver context */
 72:   KSP                ksp;
 73:   AppCtx             user;            /* user-defined work context */
 74:   PetscReal          zero = 0.0;      /* lower bound on all variables */

 76:   /* Initialize PETSC and TAO */
 77:   PetscInitialize(&argc, &argv,(char *)0,help);

 79:   /* Set the default values for the problem parameters */
 80:   user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
 81:   testgetdiag = PETSC_FALSE;

 83:   /* Check for any command line arguments that override defaults */
 84:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.nx,&flg);
 85:   PetscOptionsGetInt(NULL,NULL,"-my",&user.ny,&flg);
 86:   PetscOptionsGetReal(NULL,NULL,"-ecc",&user.ecc,&flg);
 87:   PetscOptionsGetReal(NULL,NULL,"-b",&user.b,&flg);
 88:   PetscOptionsGetBool(NULL,NULL,"-test_getdiagonal",&testgetdiag,NULL);

 90:   PetscPrintf(PETSC_COMM_WORLD,"\n---- Journal Bearing Problem SHB-----\n");
 91:   PetscPrintf(PETSC_COMM_WORLD,"mx: %D,  my: %D,  ecc: %g \n\n",user.nx,user.ny,(double)user.ecc);

 93:   /* Let Petsc determine the grid division */
 94:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 96:   /*
 97:      A two dimensional distributed array will help define this problem,
 98:      which derives from an elliptic PDE on two dimensional domain.  From
 99:      the distributed array, Create the vectors.
100:   */
101:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.nx,user.ny,Nx,Ny,1,1,NULL,NULL,&user.dm);
102:   DMSetFromOptions(user.dm);
103:   DMSetUp(user.dm);

105:   /*
106:      Extract global and local vectors from DM; the vector user.B is
107:      used solely as work space for the evaluation of the function,
108:      gradient, and Hessian.  Duplicate for remaining vectors that are
109:      the same types.
110:   */
111:   DMCreateGlobalVector(user.dm,&x); /* Solution */
112:   VecDuplicate(x,&user.B); /* Linear objective */

114:   /*  Create matrix user.A to store quadratic, Create a local ordering scheme. */
115:   VecGetLocalSize(x,&m);
116:   DMCreateMatrix(user.dm,&user.A);

118:   if (testgetdiag) {
119:     MatSetOperation(user.A,MATOP_GET_DIAGONAL,NULL);
120:   }

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

125:   /* The TAO code begins here */

127:   /*
128:      Create the optimization solver
129:      Suitable methods: TAOGPCG, TAOBQPIP, TAOTRON, TAOBLMVM
130:   */
131:   TaoCreate(PETSC_COMM_WORLD,&tao);
132:   TaoSetType(tao,TAOBLMVM);

134:   /* Set the initial vector */
135:   VecSet(x, zero);
136:   TaoSetSolution(tao,x);

138:   /* Set the user function, gradient, hessian evaluation routines and data structures */
139:   TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void*) &user);

141:   TaoSetHessian(tao,user.A,user.A,FormHessian,(void*)&user);

143:   /* Set a routine that defines the bounds */
144:   VecDuplicate(x,&xl);
145:   VecDuplicate(x,&xu);
146:   VecSet(xl, zero);
147:   VecSet(xu, d1000);
148:   TaoSetVariableBounds(tao,xl,xu);

150:   TaoGetKSP(tao,&ksp);
151:   if (ksp) {
152:     KSPSetType(ksp,KSPCG);
153:   }

155:   PetscOptionsHasName(NULL,NULL,"-testmonitor",&flg);
156:   if (flg) {
157:     TaoSetMonitor(tao,Monitor,&user,NULL);
158:   }
159:   PetscOptionsHasName(NULL,NULL,"-testconvergence",&flg);
160:   if (flg) {
161:     TaoSetConvergenceTest(tao,ConvergenceTest,&user);
162:   }

164:   /* Check for any tao command line options */
165:   TaoSetFromOptions(tao);

167:   /* Solve the bound constrained problem */
168:   TaoSolve(tao);

170:   /* Free PETSc data structures */
171:   VecDestroy(&x);
172:   VecDestroy(&xl);
173:   VecDestroy(&xu);
174:   MatDestroy(&user.A);
175:   VecDestroy(&user.B);

177:   /* Free TAO data structures */
178:   TaoDestroy(&tao);
179:   DMDestroy(&user.dm);
180:   PetscFinalize();
181:   return 0;
182: }

184: static PetscReal p(PetscReal xi, PetscReal ecc)
185: {
186:   PetscReal t=1.0+ecc*PetscCosScalar(xi);
187:   return (t*t*t);
188: }

190: PetscErrorCode ComputeB(AppCtx* user)
191: {
192:   PetscInt       i,j,k;
193:   PetscInt       nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
194:   PetscReal      two=2.0, pi=4.0*atan(1.0);
195:   PetscReal      hx,hy,ehxhy;
196:   PetscReal      temp,*b;
197:   PetscReal      ecc=user->ecc;

199:   nx=user->nx;
200:   ny=user->ny;
201:   hx=two*pi/(nx+1.0);
202:   hy=two*user->b/(ny+1.0);
203:   ehxhy = ecc*hx*hy;

205:   /*
206:      Get local grid boundaries
207:   */
208:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
209:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

211:   /* Compute the linear term in the objective function */
212:   VecGetArray(user->B,&b);
213:   for (i=xs; i<xs+xm; i++) {
214:     temp=PetscSinScalar((i+1)*hx);
215:     for (j=ys; j<ys+ym; j++) {
216:       k=xm*(j-ys)+(i-xs);
217:       b[k]=  - ehxhy*temp;
218:     }
219:   }
220:   VecRestoreArray(user->B,&b);
221:   PetscLogFlops(5.0*xm*ym+3.0*xm);
222:   return 0;
223: }

225: PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *fcn,Vec G,void *ptr)
226: {
227:   AppCtx*        user=(AppCtx*)ptr;
228:   PetscInt       i,j,k,kk;
229:   PetscInt       col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
230:   PetscReal      one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
231:   PetscReal      hx,hy,hxhy,hxhx,hyhy;
232:   PetscReal      xi,v[5];
233:   PetscReal      ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
234:   PetscReal      vmiddle, vup, vdown, vleft, vright;
235:   PetscReal      tt,f1,f2;
236:   PetscReal      *x,*g,zero=0.0;
237:   Vec            localX;

239:   nx=user->nx;
240:   ny=user->ny;
241:   hx=two*pi/(nx+1.0);
242:   hy=two*user->b/(ny+1.0);
243:   hxhy=hx*hy;
244:   hxhx=one/(hx*hx);
245:   hyhy=one/(hy*hy);

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

249:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
250:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

252:   VecSet(G, zero);
253:   /*
254:     Get local grid boundaries
255:   */
256:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
257:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

259:   VecGetArray(localX,&x);
260:   VecGetArray(G,&g);

262:   for (i=xs; i< xs+xm; i++) {
263:     xi=(i+1)*hx;
264:     trule1=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
265:     trule2=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
266:     trule3=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
267:     trule4=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
268:     trule5=trule1; /* L(i,j-1) */
269:     trule6=trule2; /* U(i,j+1) */

271:     vdown=-(trule5+trule2)*hyhy;
272:     vleft=-hxhx*(trule2+trule4);
273:     vright= -hxhx*(trule1+trule3);
274:     vup=-hyhy*(trule1+trule6);
275:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);

277:     for (j=ys; j<ys+ym; j++) {

279:       row=(j-gys)*gxm + (i-gxs);
280:        v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;

282:        k=0;
283:        if (j>gys) {
284:          v[k]=vdown; col[k]=row - gxm; k++;
285:        }

287:        if (i>gxs) {
288:          v[k]= vleft; col[k]=row - 1; k++;
289:        }

291:        v[k]= vmiddle; col[k]=row; k++;

293:        if (i+1 < gxs+gxm) {
294:          v[k]= vright; col[k]=row+1; k++;
295:        }

297:        if (j+1 <gys+gym) {
298:          v[k]= vup; col[k] = row+gxm; k++;
299:        }
300:        tt=0;
301:        for (kk=0;kk<k;kk++) {
302:          tt+=v[kk]*x[col[kk]];
303:        }
304:        row=(j-ys)*xm + (i-xs);
305:        g[row]=tt;

307:      }

309:   }

311:   VecRestoreArray(localX,&x);
312:   VecRestoreArray(G,&g);

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

316:   VecDot(X,G,&f1);
317:   VecDot(user->B,X,&f2);
318:   VecAXPY(G, one, user->B);
319:   *fcn = f1/2.0 + f2;

321:   PetscLogFlops((91 + 10.0*ym) * xm);
322:   return 0;

324: }

326: /*
327:    FormHessian computes the quadratic term in the quadratic objective function
328:    Notice that the objective function in this problem is quadratic (therefore a constant
329:    hessian).  If using a nonquadratic solver, then you might want to reconsider this function
330: */
331: PetscErrorCode FormHessian(Tao tao,Vec X,Mat hes, Mat Hpre, void *ptr)
332: {
333:   AppCtx*        user=(AppCtx*)ptr;
334:   PetscInt       i,j,k;
335:   PetscInt       col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
336:   PetscReal      one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
337:   PetscReal      hx,hy,hxhy,hxhx,hyhy;
338:   PetscReal      xi,v[5];
339:   PetscReal      ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
340:   PetscReal      vmiddle, vup, vdown, vleft, vright;
341:   PetscBool      assembled;

343:   nx=user->nx;
344:   ny=user->ny;
345:   hx=two*pi/(nx+1.0);
346:   hy=two*user->b/(ny+1.0);
347:   hxhy=hx*hy;
348:   hxhx=one/(hx*hx);
349:   hyhy=one/(hy*hy);

351:   /*
352:     Get local grid boundaries
353:   */
354:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
355:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
356:   MatAssembled(hes,&assembled);
357:   if (assembled) MatZeroEntries(hes);

359:   for (i=xs; i< xs+xm; i++) {
360:     xi=(i+1)*hx;
361:     trule1=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc)) / six; /* L(i,j) */
362:     trule2=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc)) / six; /* U(i,j) */
363:     trule3=hxhy*(p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc)) / six; /* U(i+1,j) */
364:     trule4=hxhy*(p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc)) / six; /* L(i-1,j) */
365:     trule5=trule1; /* L(i,j-1) */
366:     trule6=trule2; /* U(i,j+1) */

368:     vdown=-(trule5+trule2)*hyhy;
369:     vleft=-hxhx*(trule2+trule4);
370:     vright= -hxhx*(trule1+trule3);
371:     vup=-hyhy*(trule1+trule6);
372:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
373:     v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;

375:     for (j=ys; j<ys+ym; j++) {
376:       row=(j-gys)*gxm + (i-gxs);

378:       k=0;
379:       if (j>gys) {
380:         v[k]=vdown; col[k]=row - gxm; k++;
381:       }

383:       if (i>gxs) {
384:         v[k]= vleft; col[k]=row - 1; k++;
385:       }

387:       v[k]= vmiddle; col[k]=row; k++;

389:       if (i+1 < gxs+gxm) {
390:         v[k]= vright; col[k]=row+1; k++;
391:       }

393:       if (j+1 <gys+gym) {
394:         v[k]= vup; col[k] = row+gxm; k++;
395:       }
396:       MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES);

398:     }

400:   }

402:   /*
403:      Assemble matrix, using the 2-step process:
404:      MatAssemblyBegin(), MatAssemblyEnd().
405:      By placing code between these two statements, computations can be
406:      done while messages are in transition.
407:   */
408:   MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);
409:   MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);

411:   /*
412:     Tell the matrix we will never add a new nonzero location to the
413:     matrix. If we do it will generate an error.
414:   */
415:   MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
416:   MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);

418:   PetscLogFlops(9.0*xm*ym+49.0*xm);
419:   return 0;
420: }

422: PetscErrorCode Monitor(Tao tao, void *ctx)
423: {
424:   PetscInt           its;
425:   PetscReal          f,gnorm,cnorm,xdiff;
426:   TaoConvergedReason reason;

428:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
429:   if (!(its%5)) {
430:     PetscPrintf(PETSC_COMM_WORLD,"iteration=%D\tf=%g\n",its,(double)f);
431:   }
432:   return 0;
433: }

435: PetscErrorCode ConvergenceTest(Tao tao, void *ctx)
436: {
437:   PetscInt           its;
438:   PetscReal          f,gnorm,cnorm,xdiff;
439:   TaoConvergedReason reason;

441:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
442:   if (its == 100) {
443:     TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS);
444:   }
445:   return 0;

447: }

449: /*TEST

451:    build:
452:       requires: !complex

454:    test:
455:       args: -tao_smonitor -mx 8 -my 12 -tao_type tron -tao_gatol 1.e-5
456:       requires: !single

458:    test:
459:       suffix: 2
460:       nsize: 2
461:       args: -tao_smonitor -mx 50 -my 50 -ecc 0.99 -tao_type gpcg -tao_gatol 1.e-5
462:       requires: !single

464:    test:
465:       suffix: 3
466:       nsize: 2
467:       args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4
468:       requires: !single

470:    test:
471:       suffix: 4
472:       nsize: 2
473:       args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4 -test_getdiagonal
474:       output_file: output/jbearing2_3.out
475:       requires: !single

477:    test:
478:       suffix: 5
479:       args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_bncg_type gd -tao_gatol 1e-4
480:       requires: !single

482:    test:
483:       suffix: 6
484:       args: -tao_smonitor -mx 8 -my 12 -tao_type bncg -tao_gatol 1e-4
485:       requires: !single

487:    test:
488:       suffix: 7
489:       args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5
490:       requires: !single

492:    test:
493:       suffix: 8
494:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5
495:       requires: !single

497:    test:
498:       suffix: 9
499:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5
500:       requires: !single

502:    test:
503:       suffix: 10
504:       args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
505:       requires: !single

507:    test:
508:       suffix: 11
509:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
510:       requires: !single

512:    test:
513:       suffix: 12
514:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_bnk_max_cg_its 3
515:       requires: !single

517:    test:
518:      suffix: 13
519:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls
520:      requires: !single

522:    test:
523:      suffix: 14
524:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type blmvm
525:      requires: !single

527:    test:
528:      suffix: 15
529:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnkls -tao_bqnk_mat_type lmvmbfgs
530:      requires: !single

532:    test:
533:      suffix: 16
534:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnktr -tao_bqnk_mat_type lmvmsr1
535:      requires: !single

537:    test:
538:      suffix: 17
539:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type scalar -tao_view
540:      requires: !single

542:    test:
543:      suffix: 18
544:      args: -tao_smonitor -mx 8 -my 12 -tao_gatol 1e-4 -tao_type bqnls -tao_bqnls_mat_lmvm_scale_type none -tao_view
545:      requires: !single

547:    test:
548:      suffix: 19
549:      args: -tao_smonitor -mx 8 -my 12 -tao_type bnls -tao_gatol 1e-5 -tao_mf_hessian
550:      requires: !single

552:    test:
553:       suffix: 20
554:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntr -tao_gatol 1e-5 -tao_mf_hessian
555:       requires: !single

557:    test:
558:       suffix: 21
559:       args: -tao_smonitor -mx 8 -my 12 -tao_type bntl -tao_gatol 1e-5 -tao_mf_hessian
560:       requires: !single
561: TEST*/