Actual source code: eptorsion2.c

tao-2.1-p0 2012-07-24
  1: /* Program usage: mpirun -np <proc> eptorsion2 [-help] [all TAO options] */

  3: /* ----------------------------------------------------------------------

  5:   Elastic-plastic torsion problem.  

  7:   The elastic plastic torsion problem arises from the determination 
  8:   of the stress field on an infinitely long cylindrical bar, which is
  9:   equivalent to the solution of the following problem:

 11:   min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
 12:   
 13:   where C is the torsion angle per unit length.

 15:   The uniprocessor version of this code is eptorsion1.c; the Fortran 
 16:   version of this code is eptorsion2f.F.

 18:   This application solves the problem without calculating hessians 
 19: ---------------------------------------------------------------------- */

 21: /*
 22:   Include "tao.h" so that we can use TAO solvers.  Note that this 
 23:   file automatically includes files for lower-level support, such as those
 24:   provided by the PETSc library:
 25:      petsc.h       - base PETSc routines   petscvec.h - vectors
 26:      petscsys.h    - sysem routines        petscmat.h - matrices
 27:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 28:      petscviewer.h - viewers               petscpc.h  - preconditioners
 29:   Include "petscdm.h" so that we can use distributed arrays (DMs) for managing
 30:   the parallel mesh.
 31: */

 33: #include "taosolver.h"
 34: #include "petscdm.h"

 36: static  char help[] = 
 37: "Demonstrates use of the TAO package to solve \n\
 38: unconstrained minimization problems in parallel.  This example is based on \n\
 39: the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.\n\
 40: The command line options are:\n\
 41:   -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction\n\
 42:   -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction\n\
 43:   -par <param>, where <param> = angle of twist per unit length\n\n";

 45: /*T 
 46:    Concepts: TAO - Solving an unconstrained minimization problem
 47:    Routines: TaoInitialize(); TaoFinalize(); 
 48:    Routines: TaoCreate(); TaoSetType();
 49:    Routines: TaoSetInitialVector(); 
 50:    Routines: TaoSetObjectiveAndGradientRoutine();
 51:    Routines: TaoSetHessianRoutine(); TaoSetFromOptions();
 52:    Routines: TaoSolve();
 53:    Routines: TaoGetTerminationReason(); TaoDestroy();
 54:    Processors: n
 55: T*/ 

 57: /* 
 58:    User-defined application context - contains data needed by the 
 59:    application-provided call-back routines, FormFunction() and
 60:    FormGradient().
 61: */
 62: typedef struct {
 63:   /* parameters */
 64:    PetscInt           mx, my;         /* global discretization in x- and y-directions */
 65:    PetscReal        param;          /* nonlinearity parameter */

 67:   /* work space */
 68:    Vec           localX;         /* local vectors */
 69:    DM            dm;             /* distributed array data structure */
 70: } AppCtx;


 73: PetscErrorCode FormInitialGuess(AppCtx*, Vec);
 74: PetscErrorCode FormFunctionGradient(TaoSolver,Vec,PetscReal*,Vec,void*);
 75: PetscErrorCode FormHessian(TaoSolver,Vec,Mat*,Mat*,MatStructure*,void*);


 80: int main(int argc, char **argv) 
 81: {
 83:     Vec x;
 84:     Mat H;
 85:     PetscInt Nx, Ny;
 86:     TaoSolver tao;
 87:     TaoSolverTerminationReason reason;
 88:     PetscBool flg;
 89:     KSP ksp; PC pc;
 90:     AppCtx user;

 92:     /* Initialize PETSc, TAO */
 93:     PetscInitialize(&argc, &argv, (char *)0, help);
 94:     TaoInitialize(&argc, &argv, (char *)0, help);

 96:     /* Specify default dimension of the problem */
 97:     user.param = 5.0; user.mx = 10; user.my = 10;
 98:     Nx = Ny = PETSC_DECIDE;

100:     /* Check for any command line arguments that override defaults */
101:     PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,&flg); 
102:     PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,&flg); 
103:     PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,&flg); 

105:     PetscPrintf(PETSC_COMM_WORLD,"\n---- Elastic-Plastic Torsion Problem -----\n");
106:     PetscPrintf(PETSC_COMM_WORLD,"mx: %D     my: %D   \n\n",user.mx,user.my);  

108:     /* Set up distributed array */
109:     DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,
110:                         DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,
111:                         user.mx,user.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,
112:                         &user.dm); 

114:     /* Create vectors */
115:     DMCreateGlobalVector(user.dm,&x); 

117:     DMCreateLocalVector(user.dm,&user.localX); 

119:     /* Create Hessian */
120:     DMCreateMatrix(user.dm,MATAIJ,&H); 
121:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE); 

123:     /* The TAO code begins here */

125:     /* Create TAO solver and set desired solution method */
126:     TaoCreate(PETSC_COMM_WORLD,&tao); 
127:     TaoSetType(tao,"tao_cg"); 

129:     /* Set initial solution guess */
130:     FormInitialGuess(&user,x); 
131:     TaoSetInitialVector(tao,x); 

133:     /* Set routine for function and gradient evaluation */
134:     TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user); 

136:     TaoSetHessianRoutine(tao,H,H,FormHessian,(void*)&user); 


139:     /* Check for any TAO command line options */
140:     TaoSetFromOptions(tao); 
141:     
142:     TaoGetKSP(tao,&ksp); 
143:     if (ksp) {
144:       KSPGetPC(ksp,&pc); 
145:       PCSetType(pc,PCNONE); 
146:     }

148:     /* SOLVE THE APPLICATION */
149:     TaoSolve(tao);  

151:     /* Get information on termination */
152:     TaoGetTerminationReason(tao,&reason); 
153:     if (reason <= 0){
154:         ierr=PetscPrintf(MPI_COMM_WORLD, "Try another method! \n");
155:          
156:     }  

158:     /* Free TAO data structures */
159:     TaoDestroy(&tao); 

161:     /* Free PETSc data structures */
162:     VecDestroy(&x); 
163:     MatDestroy(&H); 

165:     VecDestroy(&user.localX); 
166:     DMDestroy(&user.dm); 


169:     /* Finalize TAO and PETSc*/
170:     TaoFinalize();
171:     PetscFinalize();
172:     
173:     return 0;
174: }


177: /* ------------------------------------------------------------------- */
180: /*
181:     FormInitialGuess - Computes an initial approximation to the solution.

183:     Input Parameters:
184: .   user - user-defined application context
185: .   X    - vector

187:     Output Parameters:
188:     X    - vector
189: */
190: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
191: {
192:   PetscErrorCode    ierr;
193:   PetscInt   i, j, k, mx = user->mx, my = user->my;
194:   PetscInt   xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
195:   PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;

198:   /* Get local mesh boundaries */
199:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
200:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

202:   /* Compute initial guess over locally owned part of mesh */
203:   xe = xs+xm;
204:   ye = ys+ym;
205:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
206:     temp = PetscMin(j+1,my-j)*hy;
207:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
208:       k   = (j-gys)*gxm + i-gxs;
209:       val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
210:       VecSetValuesLocal(X,1,&k,&val,ADD_VALUES); 
211:     }
212:   }
213:   VecAssemblyBegin(X); 
214:   VecAssemblyEnd(X); 
215:   return(0);
216: }


219: /* ------------------------------------------------------------------ */
222: /* 
223:    FormFunctionGradient - Evaluates the function and corresponding gradient.
224:     
225:    Input Parameters:
226:    tao - the TaoSolver context
227:    X   - the input vector 
228:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

230:    Output Parameters:
231:    f   - the newly evaluated function
232:    G   - the newly evaluated gradient
233: */
234: PetscErrorCode FormFunctionGradient(TaoSolver tao,Vec X,PetscReal *f,Vec G,void *ptr){

236:   AppCtx *user = (AppCtx *)ptr;
237:   PetscErrorCode    ierr;
238:   PetscInt i,j,k,ind;
239:   PetscInt xe,ye,xsm,ysm,xep,yep;
240:   PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys;
241:   PetscInt mx = user->mx, my = user->my;
242:   PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
243:   PetscReal p5 = 0.5, area, val, flin, fquad;
244:   PetscReal v,vb,vl,vr,vt,dvdx,dvdy;
245:   PetscReal hx = 1.0/(user->mx + 1); 
246:   PetscReal hy = 1.0/(user->my + 1);  
247:   Vec    localX = user->localX;


251:   /* Initialize */
252:   flin = fquad = zero;

254:   VecSet(G, zero); 
255:   /*
256:      Scatter ghost points to local vector,using the 2-step process
257:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
258:      By placing code between these two statements, computations can be
259:      done while messages are in transition.
260:   */
261:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
262:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

264:   /* Get pointer to vector data */
265:   VecGetArray(localX,&x); 

267:   /* Get local mesh boundaries */
268:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
269:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

271:   /* Set local loop dimensions */
272:   xe = xs+xm;
273:   ye = ys+ym;
274:   if (xs == 0)  xsm = xs-1;
275:   else          xsm = xs;
276:   if (ys == 0)  ysm = ys-1;
277:   else          ysm = ys;
278:   if (xe == mx) xep = xe+1;
279:   else          xep = xe;
280:   if (ye == my) yep = ye+1;
281:   else          yep = ye;

283:   /* Compute local gradient contributions over the lower triangular elements */
284:   for (j=ysm; j<ye; j++) {  /*  for (j=-1; j<my; j++) */
285:     for (i=xsm; i<xe; i++) {  /*   for (i=-1; i<mx; i++) */
286:       k = (j-gys)*gxm + i-gxs;
287:       v = zero;
288:       vr = zero;
289:       vt = zero;
290:       if (i >= 0 && j >= 0) v = x[k];
291:       if (i < mx-1 && j > -1) vr = x[k+1];
292:       if (i > -1 && j < my-1) vt = x[k+gxm];
293:       dvdx = (vr-v)/hx;
294:       dvdy = (vt-v)/hy;
295:       if (i != -1 && j != -1) {
296:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
297:         VecSetValuesLocal(G,1,&k,&val,ADD_VALUES); 
298:       }
299:       if (i != mx-1 && j != -1) {
300:         ind = k+1; val =  dvdx/hx - cdiv3;
301:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); 
302:       }
303:       if (i != -1 && j != my-1) {
304:         ind = k+gxm; val = dvdy/hy - cdiv3;
305:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); 
306:       }
307:       fquad += dvdx*dvdx + dvdy*dvdy;
308:       flin -= cdiv3 * (v + vr + vt);
309:     }
310:   }

312:   /* Compute local gradient contributions over the upper triangular elements */
313:   for (j=ys; j<yep; j++) { /*  for (j=0; j<=my; j++) */
314:     for (i=xs; i<xep; i++) {  /*   for (i=0; i<=mx; i++) */
315:       k = (j-gys)*gxm + i-gxs;
316:       vb = zero;
317:       vl = zero;
318:       v  = zero;
319:       if (i < mx && j > 0) vb = x[k-gxm];
320:       if (i > 0 && j < my) vl = x[k-1];
321:       if (i < mx && j < my) v = x[k];
322:       dvdx = (v-vl)/hx;
323:       dvdy = (v-vb)/hy;
324:       if (i != mx && j != 0) {
325:         ind = k-gxm; val = - dvdy/hy - cdiv3;
326:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); 
327:       }
328:       if (i != 0 && j != my) {
329:         ind = k-1; val =  - dvdx/hx - cdiv3;
330:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); 
331:       }
332:       if (i != mx && j != my) {
333:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
334:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES); 
335:       }
336:       fquad += dvdx*dvdx + dvdy*dvdy;
337:       flin -= cdiv3 * (vb + vl + v);
338:     }
339:   }


342:   /* Restore vector */
343:   VecRestoreArray(localX,&x); 

345:   /* Assemble gradient vector */
346:   VecAssemblyBegin(G); 
347:   VecAssemblyEnd(G); 

349:   /* Scale the gradient */
350:   area = p5*hx*hy;
351:   floc = area * (p5 * fquad + flin);
352:   VecScale(G, area); 

354:   /* Sum function contributions from all processes */
355:   (PetscErrorCode)MPI_Allreduce((void*)&floc,(void*)f,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD); 

357:   ierr=PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16); 

359:   return(0);
360: }



366: PetscErrorCode FormHessian(TaoSolver tao, Vec X, Mat *H, Mat *Hpre, MatStructure *flag, void*ctx){

368:   AppCtx *user= (AppCtx*) ctx;
370:   PetscInt i,j,k;
371:   PetscInt col[5],row;
372:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
373:   PetscReal v[5];
374:   PetscReal hx=1.0/(user->mx+1), hy=1.0/(user->my+1), hxhx=1.0/(hx*hx), hyhy=1.0/(hy*hy), area=0.5*hx*hy;
375:   Mat A=*H;

377:   /* Compute the quadratic term in the objective function */  

379:   /*
380:      Get local grid boundaries
381:   */

384:   DMDAGetCorners(user->dm,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL); 
385:   DMDAGetGhostCorners(user->dm,&gxs,&gys,PETSC_NULL,&gxm,&gym,PETSC_NULL); 

387:   for (j=ys; j<ys+ym; j++){
388:     
389:     for (i=xs; i< xs+xm; i++){

391:       row=(j-gys)*gxm + (i-gxs);

393:       k=0;
394:       if (j>gys){ 
395:         v[k]=-2*hyhy; col[k]=row - gxm; k++;
396:       }

398:       if (i>gxs){
399:         v[k]= -2*hxhx; col[k]=row - 1; k++;
400:       }

402:       v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++;

404:       if (i+1 < gxs+gxm){
405:         v[k]= -2.0*hxhx; col[k]=row+1; k++;
406:       }

408:       if (j+1 <gys+gym){
409:         v[k]= -2*hyhy; col[k] = row+gxm; k++;
410:       }

412:       MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES); 

414:     }
415:   }
416:   /* 
417:      Assemble matrix, using the 2-step process:
418:      MatAssemblyBegin(), MatAssemblyEnd().
419:      By placing code between these two statements, computations can be
420:      done while messages are in transition.
421:   */
422:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); 
423:   MatAssemblyEnd(A,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:   MatScale(A,area); 
429:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); 
430:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE); 

432:   PetscLogFlops(9*xm*ym+49*xm); 

434:   return(0);
435: }