Actual source code: jbearing2.c
petsc-3.9.4 2018-09-11
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: MatSetOperation(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
495:
496: test:
497: suffix: 5
498: args: -tao_smonitor -mx 8 -my 12 -tao_type pgd -tao_gttol 1.e-3 -tao_gatol 1e-4
499: requires: !single
500:
501: test:
502: suffix: 6
503: args: -tao_monitor -tao_view -mx 8 -my 12 -tao_type bncg -tao_gatol 1e-4
504: requires: !single
506: TEST*/