Actual source code: eptorsion2.c

  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:     AppCtx user;

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

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

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

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

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

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

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

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

122:     /* The TAO code begins here */

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

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

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

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


138:     /* Check for any TAO command line options */
139:     TaoSetFromOptions(tao); 

141:     /* SOLVE THE APPLICATION */
142:     TaoSolve(tao);  

144:     /* Get information on termination */
145:     TaoGetTerminationReason(tao,&reason); 
146:     if (reason <= 0){
147:         ierr=PetscPrintf(MPI_COMM_WORLD, "Try another method! \n");
148:          
149:     }  

151:     /* Free TAO data structures */
152:     TaoDestroy(&tao); 

154:     /* Free PETSc data structures */
155:     VecDestroy(&x); 
156:     MatDestroy(&H); 

158:     VecDestroy(&user.localX); 
159:     DMDestroy(&user.dm); 


162:     /* Finalize TAO and PETSc*/
163:     TaoFinalize();
164:     PetscFinalize();
165:     
166:     return 0;
167: }


170: /* ------------------------------------------------------------------- */
173: /*
174:     FormInitialGuess - Computes an initial approximation to the solution.

176:     Input Parameters:
177: .   user - user-defined application context
178: .   X    - vector

180:     Output Parameters:
181:     X    - vector
182: */
183: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
184: {
185:   PetscErrorCode    ierr;
186:   PetscInt   i, j, k, mx = user->mx, my = user->my;
187:   PetscInt   xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
188:   PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;

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

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


212: /* ------------------------------------------------------------------ */
215: /* 
216:    FormFunctionGradient - Evaluates the function and corresponding gradient.
217:     
218:    Input Parameters:
219:    tao - the TaoSolver context
220:    X   - the input vector 
221:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

223:    Output Parameters:
224:    f   - the newly evaluated function
225:    G   - the newly evaluated gradient
226: */
227: PetscErrorCode FormFunctionGradient(TaoSolver tao,Vec X,PetscReal *f,Vec G,void *ptr){

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


244:   /* Initialize */
245:   flin = fquad = zero;

247:   VecSet(G, zero); 
248:   /*
249:      Scatter ghost points to local vector,using the 2-step process
250:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
251:      By placing code between these two statements, computations can be
252:      done while messages are in transition.
253:   */
254:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX); 
255:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX); 

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

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

264:   /* Set local loop dimensions */
265:   xe = xs+xm;
266:   ye = ys+ym;
267:   if (xs == 0)  xsm = xs-1;
268:   else          xsm = xs;
269:   if (ys == 0)  ysm = ys-1;
270:   else          ysm = ys;
271:   if (xe == mx) xep = xe+1;
272:   else          xep = xe;
273:   if (ye == my) yep = ye+1;
274:   else          yep = ye;

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

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


335:   /* Restore vector */
336:   VecRestoreArray(localX,&x); 

338:   /* Assemble gradient vector */
339:   VecAssemblyBegin(G); 
340:   VecAssemblyEnd(G); 

342:   /* Scale the gradient */
343:   area = p5*hx*hy;
344:   floc = area * (p5 * fquad + flin);
345:   VecScale(G, area); 

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

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

352:   return(0);
353: }



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

361:   AppCtx *user= (AppCtx*) ctx;
363:   PetscInt i,j,k;
364:   PetscInt col[5],row;
365:   PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym;
366:   PetscReal v[5];
367:   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;
368:   Mat A=*H;

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

372:   /*
373:      Get local grid boundaries
374:   */

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

380:   for (j=ys; j<ys+ym; j++){
381:     
382:     for (i=xs; i< xs+xm; i++){

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

386:       k=0;
387:       if (j>gys){ 
388:         v[k]=-2*hyhy; col[k]=row - gxm; k++;
389:       }

391:       if (i>gxs){
392:         v[k]= -2*hxhx; col[k]=row - 1; k++;
393:       }

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

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

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

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

407:     }
408:   }
409:   /* 
410:      Assemble matrix, using the 2-step process:
411:      MatAssemblyBegin(), MatAssemblyEnd().
412:      By placing code between these two statements, computations can be
413:      done while messages are in transition.
414:   */
415:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); 
416:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); 
417:   /*
418:     Tell the matrix we will never add a new nonzero location to the
419:     matrix. If we do it will generate an error.
420:   */
421:   MatScale(A,area); 
422:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE); 
423:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE); 

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

427:   return(0);
428: }