Actual source code: eptorsion2.c

  1: /* Program usage: mpiexec -n <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)}

 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 "petsctao.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    - system routines        petscmat.h - matrices
 27:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 28:      petscviewer.h - viewers               petscpc.h  - preconditioners
 29:   Include "petscdmda.h" so that we can use distributed arrays (DMs) for managing
 30:   the parallel mesh.
 31: */

 33: #include <petsctao.h>
 34: #include <petscdmda.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: TaoCreate(); TaoSetType();
 48:    Routines: TaoSetSolution();
 49:    Routines: TaoSetObjectiveAndGradient();
 50:    Routines: TaoSetHessian(); TaoSetFromOptions();
 51:    Routines: TaoSolve();
 52:    Routines: TaoDestroy();
 53:    Processors: n
 54: T*/

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

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

 71: PetscErrorCode FormInitialGuess(AppCtx*, Vec);
 72: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
 73: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);

 75: int main(int argc, char **argv)
 76: {
 77:     Vec                x;
 78:     Mat                H;
 79:     PetscInt           Nx, Ny;
 80:     Tao                tao;
 81:     PetscBool          flg;
 82:     KSP                ksp;
 83:     PC                 pc;
 84:     AppCtx             user;

 86:     PetscInitialize(&argc, &argv, (char *)0, help);

 88:     /* Specify default dimension of the problem */
 89:     user.param = 5.0; user.mx = 10; user.my = 10;
 90:     Nx = Ny = PETSC_DECIDE;

 92:     /* Check for any command line arguments that override defaults */
 93:     PetscOptionsGetReal(NULL,NULL,"-par",&user.param,&flg);
 94:     PetscOptionsGetInt(NULL,NULL,"-my",&user.my,&flg);
 95:     PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,&flg);

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

100:     /* Set up distributed array */
101:     DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,NULL,NULL,&user.dm);
102:     DMSetFromOptions(user.dm);
103:     DMSetUp(user.dm);

105:     /* Create vectors */
106:     DMCreateGlobalVector(user.dm,&x);

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

110:     /* Create Hessian */
111:     DMCreateMatrix(user.dm,&H);
112:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);

114:     /* The TAO code begins here */

116:     /* Create TAO solver and set desired solution method */
117:     TaoCreate(PETSC_COMM_WORLD,&tao);
118:     TaoSetType(tao,TAOCG);

120:     /* Set initial solution guess */
121:     FormInitialGuess(&user,x);
122:     TaoSetSolution(tao,x);

124:     /* Set routine for function and gradient evaluation */
125:     TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);

127:     TaoSetHessian(tao,H,H,FormHessian,(void*)&user);

129:     /* Check for any TAO command line options */
130:     TaoSetFromOptions(tao);

132:     TaoGetKSP(tao,&ksp);
133:     if (ksp) {
134:       KSPGetPC(ksp,&pc);
135:       PCSetType(pc,PCNONE);
136:     }

138:     /* SOLVE THE APPLICATION */
139:     TaoSolve(tao);

141:     /* Free TAO data structures */
142:     TaoDestroy(&tao);

144:     /* Free PETSc data structures */
145:     VecDestroy(&x);
146:     MatDestroy(&H);

148:     VecDestroy(&user.localX);
149:     DMDestroy(&user.dm);

151:     PetscFinalize();
152:     return 0;
153: }

155: /* ------------------------------------------------------------------- */
156: /*
157:     FormInitialGuess - Computes an initial approximation to the solution.

159:     Input Parameters:
160: .   user - user-defined application context
161: .   X    - vector

163:     Output Parameters:
164:     X    - vector
165: */
166: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
167: {
168:   PetscInt       i, j, k, mx = user->mx, my = user->my;
169:   PetscInt       xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
170:   PetscReal      hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val;

172:   /* Get local mesh boundaries */
173:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
174:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

176:   /* Compute initial guess over locally owned part of mesh */
177:   xe = xs+xm;
178:   ye = ys+ym;
179:   for (j=ys; j<ye; j++) {  /*  for (j=0; j<my; j++) */
180:     temp = PetscMin(j+1,my-j)*hy;
181:     for (i=xs; i<xe; i++) {  /*  for (i=0; i<mx; i++) */
182:       k   = (j-gys)*gxm + i-gxs;
183:       val = PetscMin((PetscMin(i+1,mx-i))*hx,temp);
184:       VecSetValuesLocal(X,1,&k,&val,ADD_VALUES);
185:     }
186:   }
187:   VecAssemblyBegin(X);
188:   VecAssemblyEnd(X);
189:   return 0;
190: }

192: /* ------------------------------------------------------------------ */
193: /*
194:    FormFunctionGradient - Evaluates the function and corresponding gradient.

196:    Input Parameters:
197:    tao - the Tao context
198:    X   - the input vector
199:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

201:    Output Parameters:
202:    f   - the newly evaluated function
203:    G   - the newly evaluated gradient
204: */
205: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr)
206: {
207:   AppCtx         *user = (AppCtx *)ptr;
208:   PetscInt       i,j,k,ind;
209:   PetscInt       xe,ye,xsm,ysm,xep,yep;
210:   PetscInt       xs, ys, xm, ym, gxm, gym, gxs, gys;
211:   PetscInt       mx = user->mx, my = user->my;
212:   PetscReal      three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
213:   PetscReal      p5 = 0.5, area, val, flin, fquad;
214:   PetscReal      v,vb,vl,vr,vt,dvdx,dvdy;
215:   PetscReal      hx = 1.0/(user->mx + 1);
216:   PetscReal      hy = 1.0/(user->my + 1);
217:   Vec            localX = user->localX;

219:   /* Initialize */
220:   flin = fquad = zero;

222:   VecSet(G, zero);
223:   /*
224:      Scatter ghost points to local vector,using the 2-step process
225:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
226:      By placing code between these two statements, computations can be
227:      done while messages are in transition.
228:   */
229:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
230:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

235:   /* Get local mesh boundaries */
236:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
237:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

239:   /* Set local loop dimensions */
240:   xe = xs+xm;
241:   ye = ys+ym;
242:   if (xs == 0)  xsm = xs-1;
243:   else          xsm = xs;
244:   if (ys == 0)  ysm = ys-1;
245:   else          ysm = ys;
246:   if (xe == mx) xep = xe+1;
247:   else          xep = xe;
248:   if (ye == my) yep = ye+1;
249:   else          yep = ye;

251:   /* Compute local gradient contributions over the lower triangular elements */
252:   for (j=ysm; j<ye; j++) {  /*  for (j=-1; j<my; j++) */
253:     for (i=xsm; i<xe; i++) {  /*   for (i=-1; i<mx; i++) */
254:       k = (j-gys)*gxm + i-gxs;
255:       v = zero;
256:       vr = zero;
257:       vt = zero;
258:       if (i >= 0 && j >= 0) v = x[k];
259:       if (i < mx-1 && j > -1) vr = x[k+1];
260:       if (i > -1 && j < my-1) vt = x[k+gxm];
261:       dvdx = (vr-v)/hx;
262:       dvdy = (vt-v)/hy;
263:       if (i != -1 && j != -1) {
264:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
265:         VecSetValuesLocal(G,1,&k,&val,ADD_VALUES);
266:       }
267:       if (i != mx-1 && j != -1) {
268:         ind = k+1; val =  dvdx/hx - cdiv3;
269:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
270:       }
271:       if (i != -1 && j != my-1) {
272:         ind = k+gxm; val = dvdy/hy - cdiv3;
273:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
274:       }
275:       fquad += dvdx*dvdx + dvdy*dvdy;
276:       flin -= cdiv3 * (v + vr + vt);
277:     }
278:   }

280:   /* Compute local gradient contributions over the upper triangular elements */
281:   for (j=ys; j<yep; j++) { /*  for (j=0; j<=my; j++) */
282:     for (i=xs; i<xep; i++) {  /*   for (i=0; i<=mx; i++) */
283:       k = (j-gys)*gxm + i-gxs;
284:       vb = zero;
285:       vl = zero;
286:       v  = zero;
287:       if (i < mx && j > 0) vb = x[k-gxm];
288:       if (i > 0 && j < my) vl = x[k-1];
289:       if (i < mx && j < my) v = x[k];
290:       dvdx = (v-vl)/hx;
291:       dvdy = (v-vb)/hy;
292:       if (i != mx && j != 0) {
293:         ind = k-gxm; val = - dvdy/hy - cdiv3;
294:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
295:       }
296:       if (i != 0 && j != my) {
297:         ind = k-1; val =  - dvdx/hx - cdiv3;
298:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
299:       }
300:       if (i != mx && j != my) {
301:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
302:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
303:       }
304:       fquad += dvdx*dvdx + dvdy*dvdy;
305:       flin -= cdiv3 * (vb + vl + v);
306:     }
307:   }

309:   /* Restore vector */
310:   VecRestoreArray(localX,&x);

312:   /* Assemble gradient vector */
313:   VecAssemblyBegin(G);
314:   VecAssemblyEnd(G);

316:   /* Scale the gradient */
317:   area = p5*hx*hy;
318:   floc = area * (p5 * fquad + flin);
319:   VecScale(G, area);

321:   /* Sum function contributions from all processes */  /* TODO: Change to PetscCallMPI() */
322:   (PetscErrorCode)MPI_Allreduce((void*)&floc,(void*)f,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD);

324:   PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16);
325:   return 0;
326: }

328: PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void*ctx)
329: {
330:   AppCtx         *user= (AppCtx*) ctx;
331:   PetscInt       i,j,k;
332:   PetscInt       col[5],row;
333:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
334:   PetscReal      v[5];
335:   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;

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

339:   /*
340:      Get local grid boundaries
341:   */

343:   DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL);
344:   DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL);

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

348:     for (i=xs; i< xs+xm; i++) {

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

352:       k=0;
353:       if (j>gys) {
354:         v[k]=-2*hyhy; col[k]=row - gxm; k++;
355:       }

357:       if (i>gxs) {
358:         v[k]= -2*hxhx; col[k]=row - 1; k++;
359:       }

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

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

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

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

373:     }
374:   }
375:   /*
376:      Assemble matrix, using the 2-step process:
377:      MatAssemblyBegin(), MatAssemblyEnd().
378:      By placing code between these two statements, computations can be
379:      done while messages are in transition.
380:   */
381:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
382:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
383:   /*
384:     Tell the matrix we will never add a new nonzero location to the
385:     matrix. If we do it will generate an error.
386:   */
387:   MatScale(A,area);
388:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
389:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
390:   PetscLogFlops(9*xm*ym+49*xm);
391:   return 0;
392: }

394: /*TEST

396:    build:
397:       requires: !complex

399:    test:
400:       suffix: 1
401:       nsize: 2
402:       args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4

404: TEST*/