Actual source code: jbearing2.c

petsc-master 2018-02-18
Report Typos and Errors
  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(); TaoSetObjectiveAndGradientRoutine();
 28:    Routines: TaoSetHessianRoutine();
 29:    Routines: TaoSetVariableBounds();
 30:    Routines: TaoSetMonitor(); TaoSetConvergenceTest();
 31:    Routines: TaoSetInitialVector();
 32:    Routines: TaoSetFromOptions();
 33:    Routines: TaoSolve();
 34:    Routines: TaoDestroy();
 35:    Processors: n
 36: 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(Tao, Vec, PetscReal *,Vec,void *);
 60: static PetscErrorCode FormHessian(Tao,Vec,Mat, Mat, void *);
 61: static PetscErrorCode ComputeB(AppCtx*);
 62: static PetscErrorCode Monitor(Tao, void*);
 63: static PetscErrorCode ConvergenceTest(Tao, void*);

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

 79:   /* Initialize PETSC and TAO */
 80:   PetscInitialize( &argc, &argv,(char *)0,help );if (ierr) return ierr;

 82:   /* Set the default values for the problem parameters */
 83:   user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
 84:   testgetdiag = PETSC_FALSE;

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

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

 96:   /* Let Petsc determine the grid division */
 97:   Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;

 99:   /*
100:      A two dimensional distributed array will help define this problem,
101:      which derives from an elliptic PDE on two dimensional domain.  From
102:      the distributed array, Create the vectors.
103:   */
104:   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);
105:   DMSetFromOptions(user.dm);
106:   DMSetUp(user.dm);

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


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

122:   if (testgetdiag) {
123:     MatShellSetOperation(user.A,MATOP_GET_DIAGONAL,NULL);
124:   }

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

129:   /* The TAO code begins here */

131:   /*
132:      Create the optimization solver
133:      Suitable methods: TAOGPCG, TAOBQPIP, TAOTRON, TAOBLMVM
134:   */
135:   TaoCreate(PETSC_COMM_WORLD,&tao);
136:   TaoSetType(tao,TAOBLMVM);


139:   /* Set the initial vector */
140:   VecSet(x, zero);
141:   TaoSetInitialVector(tao,x);

143:   /* Set the user function, gradient, hessian evaluation routines and data structures */
144:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*) &user);

146:   TaoSetHessianRoutine(tao,user.A,user.A,FormHessian,(void*)&user);

148:   /* Set a routine that defines the bounds */
149:   VecDuplicate(x,&xl);
150:   VecDuplicate(x,&xu);
151:   VecSet(xl, zero);
152:   VecSet(xu, d1000);
153:   TaoSetVariableBounds(tao,xl,xu);

155:   TaoGetKSP(tao,&ksp);
156:   if (ksp) {
157:     KSPSetType(ksp,KSPCG);
158:   }

160:   PetscOptionsHasName(NULL,NULL,"-testmonitor",&flg);
161:   if (flg) {
162:     TaoSetMonitor(tao,Monitor,&user,NULL);
163:   }
164:   PetscOptionsHasName(NULL,NULL,"-testconvergence",&flg);
165:   if (flg) {
166:     TaoSetConvergenceTest(tao,ConvergenceTest,&user);
167:   }

169:   /* Check for any tao command line options */
170:   TaoSetFromOptions(tao);

172:   /* Solve the bound constrained problem */
173:   TaoSolve(tao);

175:   /* Free PETSc data structures */
176:   VecDestroy(&x);
177:   VecDestroy(&xl);
178:   VecDestroy(&xu);
179:   MatDestroy(&user.A);
180:   VecDestroy(&user.B);

182:   /* Free TAO data structures */
183:   TaoDestroy(&tao);
184:   DMDestroy(&user.dm);
185:   PetscFinalize();
186:   return ierr;
187: }


190: static PetscReal p(PetscReal xi, PetscReal ecc)
191: {
192:   PetscReal t=1.0+ecc*PetscCosScalar(xi);
193:   return (t*t*t);
194: }

196: PetscErrorCode ComputeB(AppCtx* user)
197: {
199:   PetscInt       i,j,k;
200:   PetscInt       nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
201:   PetscReal      two=2.0, pi=4.0*atan(1.0);
202:   PetscReal      hx,hy,ehxhy;
203:   PetscReal      temp,*b;
204:   PetscReal      ecc=user->ecc;

206:   nx=user->nx;
207:   ny=user->ny;
208:   hx=two*pi/(nx+1.0);
209:   hy=two*user->b/(ny+1.0);
210:   ehxhy = ecc*hx*hy;


213:   /*
214:      Get local grid boundaries
215:   */
216:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
217:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

219:   /* Compute the linear term in the objective function */
220:   VecGetArray(user->B,&b);
221:   for (i=xs; i<xs+xm; i++){
222:     temp=PetscSinScalar((i+1)*hx);
223:     for (j=ys; j<ys+ym; j++){
224:       k=xm*(j-ys)+(i-xs);
225:       b[k]=  - ehxhy*temp;
226:     }
227:   }
228:   VecRestoreArray(user->B,&b);
229:   PetscLogFlops(5*xm*ym+3*xm);

231:   return 0;
232: }

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

249:   nx=user->nx;
250:   ny=user->ny;
251:   hx=two*pi/(nx+1.0);
252:   hy=two*user->b/(ny+1.0);
253:   hxhy=hx*hy;
254:   hxhx=one/(hx*hx);
255:   hyhy=one/(hy*hy);

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

259:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
260:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

262:   VecSet(G, zero);
263:   /*
264:     Get local grid boundaries
265:   */
266:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
267:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

269:   VecGetArray(localX,&x);
270:   VecGetArray(G,&g);

272:   for (i=xs; i< xs+xm; i++){
273:     xi=(i+1)*hx;
274:     trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
275:     trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
276:     trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
277:     trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
278:     trule5=trule1; /* L(i,j-1) */
279:     trule6=trule2; /* U(i,j+1) */

281:     vdown=-(trule5+trule2)*hyhy;
282:     vleft=-hxhx*(trule2+trule4);
283:     vright= -hxhx*(trule1+trule3);
284:     vup=-hyhy*(trule1+trule6);
285:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);

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

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

292:        k=0;
293:        if (j>gys){
294:          v[k]=vdown; col[k]=row - gxm; k++;
295:        }

297:        if (i>gxs){
298:          v[k]= vleft; col[k]=row - 1; k++;
299:        }

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

303:        if (i+1 < gxs+gxm){
304:          v[k]= vright; col[k]=row+1; k++;
305:        }

307:        if (j+1 <gys+gym){
308:          v[k]= vup; col[k] = row+gxm; k++;
309:        }
310:        tt=0;
311:        for (kk=0;kk<k;kk++){
312:          tt+=v[kk]*x[col[kk]];
313:        }
314:        row=(j-ys)*xm + (i-xs);
315:        g[row]=tt;

317:      }

319:   }

321:   VecRestoreArray(localX,&x);
322:   VecRestoreArray(G,&g);

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

326:   VecDot(X,G,&f1);
327:   VecDot(user->B,X,&f2);
328:   VecAXPY(G, one, user->B);
329:   *fcn = f1/2.0 + f2;


332:   PetscLogFlops((91 + 10*ym) * xm);
333:   return 0;

335: }


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

356:   nx=user->nx;
357:   ny=user->ny;
358:   hx=two*pi/(nx+1.0);
359:   hy=two*user->b/(ny+1.0);
360:   hxhy=hx*hy;
361:   hxhx=one/(hx*hx);
362:   hyhy=one/(hy*hy);

364:   /*
365:     Get local grid boundaries
366:   */
367:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
368:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);
369:   MatAssembled(hes,&assembled);
370:   if (assembled){MatZeroEntries(hes);}

372:   for (i=xs; i< xs+xm; i++){
373:     xi=(i+1)*hx;
374:     trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
375:     trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
376:     trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
377:     trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
378:     trule5=trule1; /* L(i,j-1) */
379:     trule6=trule2; /* U(i,j+1) */

381:     vdown=-(trule5+trule2)*hyhy;
382:     vleft=-hxhx*(trule2+trule4);
383:     vright= -hxhx*(trule1+trule3);
384:     vup=-hyhy*(trule1+trule6);
385:     vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
386:     v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;

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

391:       k=0;
392:       if (j>gys){
393:         v[k]=vdown; col[k]=row - gxm; k++;
394:       }

396:       if (i>gxs){
397:         v[k]= vleft; col[k]=row - 1; k++;
398:       }

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

402:       if (i+1 < gxs+gxm){
403:         v[k]= vright; col[k]=row+1; k++;
404:       }

406:       if (j+1 <gys+gym){
407:         v[k]= vup; col[k] = row+gxm; k++;
408:       }
409:       MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES);

411:     }

413:   }

415:   /*
416:      Assemble matrix, using the 2-step process:
417:      MatAssemblyBegin(), MatAssemblyEnd().
418:      By placing code between these two statements, computations can be
419:      done while messages are in transition.
420:   */
421:   MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY);
422:   MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY);

424:   /*
425:     Tell the matrix we will never add a new nonzero location to the
426:     matrix. If we do it will generate an error.
427:   */
428:   MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
429:   MatSetOption(hes,MAT_SYMMETRIC,PETSC_TRUE);

431:   PetscLogFlops(9*xm*ym+49*xm);
432:   MatNorm(hes,NORM_1,&hx);
433:   return 0;
434: }

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

444:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
445:   if (!(its%5)) {
446:     PetscPrintf(PETSC_COMM_WORLD,"iteration=%D\tf=%g\n",its,(double)f);
447:   }
448:   return(0);
449: }

451: PetscErrorCode ConvergenceTest(Tao tao, void *ctx)
452: {
453:   PetscErrorCode     ierr;
454:   PetscInt           its;
455:   PetscReal          f,gnorm,cnorm,xdiff;
456:   TaoConvergedReason reason;

459:   TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
460:   if (its == 100) {
461:     TaoSetConvergedReason(tao,TAO_DIVERGED_MAXITS);
462:   }
463:   return(0);

465: }


468: /*TEST

470:    build:
471:       requires: !complex

473:    test:
474:       args: -tao_smonitor -mx 8 -my 12 -tao_type tron -tao_gttol 1.e-5
475:       requires: !single

477:    test:
478:       suffix: 2
479:       nsize: 2
480:       args: -tao_smonitor -mx 50 -my 50 -ecc 0.99 -tao_type gpcg -tao_gttol 1.e-5
481:       requires: !single

483:    test:
484:       suffix: 3
485:       nsize: 2
486:       args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4
487:       requires: !single

489:    test:
490:       suffix: 4
491:       nsize: 2
492:       args: -tao_smonitor -mx 10 -my 16 -ecc 0.9 -tao_type bqpip -tao_gatol 1.e-4 -test_getdiagonal
493:       output_file: output/jbearing2_3.out
494:       requires: !single

496: TEST*/