Actual source code: jbearing2.c
1: /*$Id: jbearing2.c 1.78 05/05/12 13:59:03-05:00 sarich@zorak.(none) $*/
3: /*
4: Include "tao.h" so we can use TAO solvers with PETSc support.
5: Include "petscda.h" so that we can use distributed arrays (DAs) for managing
6: the parallel mesh.
7: */
9: #include "petscda.h"
10: #include tao.h
11: #include <math.h> /* for cos() sin(0), and atan() */
13: static char help[]=
14: "This example demonstrates use of the TAO package to \n\
15: solve a bound constrained minimization problem. This example is based on \n\
16: the problem DPJB from the MINPACK-2 test suite. This pressure journal \n\
17: bearing problem is an example of elliptic variational problem defined over \n\
18: a two dimensional rectangle. By discretizing the domain into triangular \n\
19: elements, the pressure surrounding the journal bearing is defined as the \n\
20: minimum of a quadratic function whose variables are bounded below by zero.\n\
21: The command line options are:\n\
22: -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
23: -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
24: \n";
26: /*T
27: Concepts: TAO - Solving a bound constrained minimization problem
28: Routines: TaoInitialize(); TaoFinalize();
29: Routines: TaoApplicationCreate(); TaoAppDestroy();
30: Routines: TaoAppSetObjectiveAndGradientRoutine(); TaoAppSetHessianRoutine();
31: Routines: TaoAppSetInitialSolutionVec(); TaoAppSetHessianMat();
32: Routines: TaoCreate(); TaoDestroy();
33: Routines: TaoSetOptions(); TaoGetKSP()
34: Routines: TaoSolveApplication(); TaoGetTerminationReason();
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: int nx,ny; /* discretization in x, y directions */
49: /* Working space */
50: DA da; /* 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 int FormFunctionGradient(TAO_APPLICATION, Vec, double *,Vec,void *);
58: static int FormHessian(TAO_APPLICATION,Vec,Mat *, Mat *, MatStructure *, void *);
59: static int ComputeB(AppCtx*);
60: static int ComputeBounds(TAO_APPLICATION, Vec, Vec, void *);
64: int main( int argc, char **argv )
65: {
66: int info; /* used to check for functions returning nonzeros */
67: int Nx, Ny; /* number of processors in x- and y- directions */
68: int m, N; /* number of local and global elements in vectors */
69: Vec x; /* variables vector */
70: PetscTruth flg; /* A return variable when checking for user options */
71: TAO_SOLVER tao; /* TAO_SOLVER solver context */
72: TaoMethod method = "tao_gpcg";/* default minimization method */
73: TAO_APPLICATION jbearingapp; /* The PETSc application */
74: ISLocalToGlobalMapping isltog; /* local-to-global mapping object */
75: int nloc; /* The number of local elements */
76: int *ltog; /* mapping of local elements to global elements */
77: TaoTerminateReason reason;
78: AppCtx user; /* user-defined work context */
79: PetscReal zero=0.0; /* lower bound on all variables */
80: KSP ksp;
82:
83: /* Initialize PETSC and TAO */
84: PetscInitialize( &argc, &argv,(char *)0,help );
85: TaoInitialize( &argc, &argv,(char *)0,help );
87: /* Set the default values for the problem parameters */
88: user.nx = 50; user.ny = 50; user.ecc = 0.1; user.b = 10.0;
90: /* Check for any command line arguments that override defaults */
91: info = PetscOptionsGetInt(TAO_NULL,"-mx",&user.nx,&flg); CHKERRQ(info);
92: info = PetscOptionsGetInt(TAO_NULL,"-my",&user.ny,&flg); CHKERRQ(info);
93: info = PetscOptionsGetReal(TAO_NULL,"-ecc",&user.ecc,&flg); CHKERRQ(info);
94: info = PetscOptionsGetReal(TAO_NULL,"-b",&user.b,&flg); CHKERRQ(info);
97: PetscPrintf(MPI_COMM_WORLD,"\n---- Journal Bearing Problem -----\n");
98: PetscPrintf(MPI_COMM_WORLD,"mx: %d, my: %d, ecc: %4.3f \n\n",
99: user.nx,user.ny,user.ecc);
101: /* Calculate any derived values from parameters */
102: N = user.nx*user.ny;
104: /* Let Petsc determine the grid division */
105: Nx = PETSC_DECIDE; Ny = PETSC_DECIDE;
107: /*
108: A two dimensional distributed array will help define this problem,
109: which derives from an elliptic PDE on two dimensional domain. From
110: the distributed array, Create the vectors.
111: */
112: info = DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.nx,
113: user.ny,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.da); CHKERRQ(info);
115: /*
116: Extract global and local vectors from DA; 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: info = DACreateGlobalVector(user.da,&x); CHKERRQ(info); /* Solution */
122: info = VecDuplicate(x,&user.B); CHKERRQ(info); /* Linear objective */
125: /* Create matrix user.A to store quadratic, Create a local ordering scheme. */
126: info = VecGetLocalSize(x,&m); CHKERRQ(info);
127: info = MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,5,TAO_NULL,3,TAO_NULL,&user.A); CHKERRQ(info);
128:
129: info = DAGetGlobalIndices(user.da,&nloc,<og); CHKERRQ(info);
130: info = ISLocalToGlobalMappingCreate(MPI_COMM_SELF,nloc,ltog,&isltog);
131: CHKERRQ(info);
132: info = MatSetLocalToGlobalMapping(user.A,isltog); CHKERRQ(info);
133: info = ISLocalToGlobalMappingDestroy(isltog); CHKERRQ(info);
136: /* User defined function -- compute linear term of quadratic */
137: info = ComputeB(&user); CHKERRQ(info);
140: /* The TAO code begins here */
142: /*
143: Create the optimization solver, Petsc application
144: Suitable methods: "tao_gpcg","tao_bqpip","tao_tron","tao_blmvm"
145: */
146: info = TaoCreate(MPI_COMM_WORLD,method,&tao); CHKERRQ(info);
147: info = TaoApplicationCreate(MPI_COMM_WORLD,&jbearingapp); CHKERRQ(info);
149: /* Set the initial vector */
150: info = VecSet(x, zero); CHKERRQ(info);
151: info = TaoAppSetInitialSolutionVec(jbearingapp,x); CHKERRQ(info);
153: /* Set the user function, gradient, hessian evaluation routines and data structures */
154: info = TaoAppSetObjectiveAndGradientRoutine(jbearingapp,FormFunctionGradient,(void*) &user);
155: CHKERRQ(info);
156:
157: info = TaoAppSetHessianMat(jbearingapp,user.A,user.A); CHKERRQ(info);
158: info = TaoAppSetHessianRoutine(jbearingapp,FormHessian,(void*)&user);
159: CHKERRQ(info);
161: /* Set a routine that defines the bounds */
162: info = TaoAppSetVariableBoundsRoutine(jbearingapp,ComputeBounds,(void*)&user); CHKERRQ(info);
164: info = TaoGetKSP(tao,&ksp); CHKERRQ(info);
165: if (ksp) { /* Modify the PETSc KSP structure */
166: info = KSPSetType(ksp,KSPCG); CHKERRQ(info);
167: }
170: /* Check for any tao command line options */
171: info = TaoSetOptions(jbearingapp,tao); CHKERRQ(info);
173: /* Solve the bound constrained problem */
174: info = TaoSolveApplication(jbearingapp,tao); CHKERRQ(info);
176: info = TaoGetTerminationReason(tao,&reason); CHKERRQ(info);
177: if (reason <= 0)
178: PetscPrintf(MPI_COMM_WORLD,"Try a different TAO method, adjust some parameters, or check the function evaluation routines\n");
180: /* Free TAO data structures */
181: info = TaoDestroy(tao); CHKERRQ(info);
182: info = TaoAppDestroy(jbearingapp); CHKERRQ(info);
184: /* Free PETSc data structures */
185: info = VecDestroy(x); CHKERRQ(info);
186: info = MatDestroy(user.A); CHKERRQ(info);
187: info = VecDestroy(user.B); CHKERRQ(info);
188: info = DADestroy(user.da); CHKERRQ(info);
190: TaoFinalize();
191: PetscFinalize();
193: return 0;
194: }
198: static int ComputeBounds(TAO_APPLICATION taoapp, Vec xl, Vec xu, void *ctx){
199: int info;
200: PetscReal zero=0.0, d1000=1000;
201: /* Set the upper and lower bounds */
202: info = VecSet(xl, zero); CHKERRQ(info);
203: info = VecSet(xu, d1000); CHKERRQ(info);
204: return 0;
205: }
207: static PetscReal p(PetscReal xi, PetscReal ecc)
208: {
209: PetscReal t=1.0+ecc*cos(xi);
210: return (t*t*t);
211: }
215: int ComputeB(AppCtx* user)
216: {
217: int i,j,k,info;
218: int nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
219: PetscReal two=2.0, pi=4.0*atan(1.0);
220: PetscReal hx,hy,ehxhy;
221: PetscReal temp,*b;
222: PetscReal ecc=user->ecc;
224: nx=user->nx;
225: ny=user->ny;
226: hx=two*pi/(nx+1.0);
227: hy=two*user->b/(ny+1.0);
228: ehxhy = ecc*hx*hy;
231: /*
232: Get local grid boundaries
233: */
234: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
235: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
236:
238: /* Compute the linear term in the objective function */
239: info = VecGetArray(user->B,&b); CHKERRQ(info);
240: for (i=xs; i<xs+xm; i++){
241: temp=sin((i+1)*hx);
242: for (j=ys; j<ys+ym; j++){
243: k=xm*(j-ys)+(i-xs);
244: b[k]= - ehxhy*temp;
245: }
246: }
247: info = VecRestoreArray(user->B,&b); CHKERRQ(info);
248: info = PetscLogFlops(5*xm*ym+3*xm); CHKERRQ(info);
250: return 0;
251: }
255: int FormFunctionGradient(TAO_APPLICATION taoapp, Vec X, double *fcn,Vec G,void *ptr)
256: {
257: AppCtx* user=(AppCtx*)ptr;
258: int i,j,k,kk,info;
259: int col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
260: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
261: PetscReal hx,hy,hxhy,hxhx,hyhy;
262: PetscReal xi,v[5];
263: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
264: PetscReal vmiddle, vup, vdown, vleft, vright;
265: PetscReal tt,f1,f2;
266: PetscReal *x,*g,zero=0.0;
267: Vec localX;
269: nx=user->nx;
270: ny=user->ny;
271: hx=two*pi/(nx+1.0);
272: hy=two*user->b/(ny+1.0);
273: hxhy=hx*hy;
274: hxhx=one/(hx*hx);
275: hyhy=one/(hy*hy);
277: info = DAGetLocalVector(user->da,&localX);CHKERRQ(info);
279: info = DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
280: info = DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX); CHKERRQ(info);
282: info = VecSet(G, zero); CHKERRQ(info);
283: /*
284: Get local grid boundaries
285: */
286: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
287: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
288:
289: info = VecGetArray(localX,&x); CHKERRQ(info);
290: info = VecGetArray(G,&g); CHKERRQ(info);
292: for (i=xs; i< xs+xm; i++){
293: xi=(i+1)*hx;
294: trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
295: trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
296: trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
297: trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
298: trule5=trule1; /* L(i,j-1) */
299: trule6=trule2; /* U(i,j+1) */
301: vdown=-(trule5+trule2)*hyhy;
302: vleft=-hxhx*(trule2+trule4);
303: vright= -hxhx*(trule1+trule3);
304: vup=-hyhy*(trule1+trule6);
305: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
307: for (j=ys; j<ys+ym; j++){
308:
309: row=(j-gys)*gxm + (i-gxs);
310: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
311:
312: k=0;
313: if (j>gys){
314: v[k]=vdown; col[k]=row - gxm; k++;
315: }
316:
317: if (i>gxs){
318: v[k]= vleft; col[k]=row - 1; k++;
319: }
321: v[k]= vmiddle; col[k]=row; k++;
322:
323: if (i+1 < gxs+gxm){
324: v[k]= vright; col[k]=row+1; k++;
325: }
326:
327: if (j+1 <gys+gym){
328: v[k]= vup; col[k] = row+gxm; k++;
329: }
330: tt=0;
331: for (kk=0;kk<k;kk++){
332: tt+=v[kk]*x[col[kk]];
333: }
334: row=(j-ys)*xm + (i-xs);
335: g[row]=tt;
337: }
339: }
341: info = VecRestoreArray(localX,&x); CHKERRQ(info);
342: info = VecRestoreArray(G,&g); CHKERRQ(info);
344: info = DARestoreLocalVector(user->da,&localX); CHKERRQ(info);
346: info = VecDot(X,G,&f1); CHKERRQ(info);
347: info = VecDot(user->B,X,&f2); CHKERRQ(info);
348: info = VecAXPY(G, one, user->B); CHKERRQ(info);
349: *fcn = f1/2.0 + f2;
351: info = PetscLogFlops((91 + 10*ym) * xm); CHKERRQ(info);
352: return 0;
354: }
360: /*
361: FormHessian computes the quadratic term in the quadratic objective function
362: Notice that the objective function in this problem is quadratic (therefore a constant
363: hessian). If using a nonquadratic solver, then you might want to reconsider this function
364: */
365: int FormHessian(TAO_APPLICATION taoapp,Vec X,Mat *H, Mat *Hpre, MatStructure *flg, void *ptr)
366: {
367: AppCtx* user=(AppCtx*)ptr;
368: int i,j,k,info;
369: int col[5],row,nx,ny,xs,xm,gxs,gxm,ys,ym,gys,gym;
370: PetscReal one=1.0, two=2.0, six=6.0,pi=4.0*atan(1.0);
371: PetscReal hx,hy,hxhy,hxhx,hyhy;
372: PetscReal xi,v[5];
373: PetscReal ecc=user->ecc, trule1,trule2,trule3,trule4,trule5,trule6;
374: PetscReal vmiddle, vup, vdown, vleft, vright;
375: Mat hes=*H;
376: PetscTruth assembled;
378: nx=user->nx;
379: ny=user->ny;
380: hx=two*pi/(nx+1.0);
381: hy=two*user->b/(ny+1.0);
382: hxhy=hx*hy;
383: hxhx=one/(hx*hx);
384: hyhy=one/(hy*hy);
386: *flg=SAME_NONZERO_PATTERN;
387: /*
388: Get local grid boundaries
389: */
390: info = DAGetCorners(user->da,&xs,&ys,TAO_NULL,&xm,&ym,TAO_NULL); CHKERRQ(info);
391: info = DAGetGhostCorners(user->da,&gxs,&gys,TAO_NULL,&gxm,&gym,TAO_NULL); CHKERRQ(info);
392:
393: info = MatAssembled(hes,&assembled); CHKERRQ(info);
394: if (assembled){info = MatZeroEntries(hes); CHKERRQ(info);}
396: for (i=xs; i< xs+xm; i++){
397: xi=(i+1)*hx;
398: trule1=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi,ecc) ) / six; /* L(i,j) */
399: trule2=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi,ecc) ) / six; /* U(i,j) */
400: trule3=hxhy*( p(xi,ecc) + p(xi+hx,ecc) + p(xi+hx,ecc) ) / six; /* U(i+1,j) */
401: trule4=hxhy*( p(xi,ecc) + p(xi-hx,ecc) + p(xi-hx,ecc) ) / six; /* L(i-1,j) */
402: trule5=trule1; /* L(i,j-1) */
403: trule6=trule2; /* U(i,j+1) */
405: vdown=-(trule5+trule2)*hyhy;
406: vleft=-hxhx*(trule2+trule4);
407: vright= -hxhx*(trule1+trule3);
408: vup=-hyhy*(trule1+trule6);
409: vmiddle=(hxhx)*(trule1+trule2+trule3+trule4)+hyhy*(trule1+trule2+trule5+trule6);
410: v[0]=0; v[1]=0; v[2]=0; v[3]=0; v[4]=0;
412: for (j=ys; j<ys+ym; j++){
413: row=(j-gys)*gxm + (i-gxs);
414:
415: k=0;
416: if (j>gys){
417: v[k]=vdown; col[k]=row - gxm; k++;
418: }
419:
420: if (i>gxs){
421: v[k]= vleft; col[k]=row - 1; k++;
422: }
424: v[k]= vmiddle; col[k]=row; k++;
425:
426: if (i+1 < gxs+gxm){
427: v[k]= vright; col[k]=row+1; k++;
428: }
429:
430: if (j+1 <gys+gym){
431: v[k]= vup; col[k] = row+gxm; k++;
432: }
433: info = MatSetValuesLocal(hes,1,&row,k,col,v,INSERT_VALUES); CHKERRQ(info);
434:
435: }
437: }
439: /*
440: Assemble matrix, using the 2-step process:
441: MatAssemblyBegin(), MatAssemblyEnd().
442: By placing code between these two statements, computations can be
443: done while messages are in transition.
444: */
445: info = MatAssemblyBegin(hes,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
446: info = MatAssemblyEnd(hes,MAT_FINAL_ASSEMBLY); CHKERRQ(info);
448: /*
449: Tell the matrix we will never add a new nonzero location to the
450: matrix. If we do it will generate an error.
451: */
452: info = MatSetOption(hes,MAT_NEW_NONZERO_LOCATION_ERR); CHKERRQ(info);
453: info = MatSetOption(hes,MAT_SYMMETRIC); CHKERRQ(info);
455: info = PetscLogFlops(9*xm*ym+49*xm); CHKERRQ(info);
457: return 0;
458: }