Actual source code: eptorsion2.c

petsc-master 2017-12-14
Report Typos and Errors
  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    - sysem 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: TaoSetInitialVector();
 49:    Routines: TaoSetObjectiveAndGradientRoutine();
 50:    Routines: TaoSetHessianRoutine(); 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;


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


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

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

 91:     /* Specify default dimension of the problem */
 92:     user.param = 5.0; user.mx = 10; user.my = 10;
 93:     Nx = Ny = PETSC_DECIDE;

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

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

103:     /* Set up distributed array */
104:     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);
105:     DMSetFromOptions(user.dm);
106:     DMSetUp(user.dm);

108:     /* Create vectors */
109:     DMCreateGlobalVector(user.dm,&x);

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

113:     /* Create Hessian */
114:     DMCreateMatrix(user.dm,&H);
115:     MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE);

117:     /* The TAO code begins here */

119:     /* Create TAO solver and set desired solution method */
120:     TaoCreate(PETSC_COMM_WORLD,&tao);
121:     TaoSetType(tao,TAOCG);

123:     /* Set initial solution guess */
124:     FormInitialGuess(&user,x);
125:     TaoSetInitialVector(tao,x);

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

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


133:     /* Check for any TAO command line options */
134:     TaoSetFromOptions(tao);

136:     TaoGetKSP(tao,&ksp);
137:     if (ksp) {
138:       KSPGetPC(ksp,&pc);
139:       PCSetType(pc,PCNONE);
140:     }

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

145:     /* Free TAO data structures */
146:     TaoDestroy(&tao);

148:     /* Free PETSc data structures */
149:     VecDestroy(&x);
150:     MatDestroy(&H);

152:     VecDestroy(&user.localX);
153:     DMDestroy(&user.dm);

155:     PetscFinalize();
156:     return 0;
157: }


160: /* ------------------------------------------------------------------- */
161: /*
162:     FormInitialGuess - Computes an initial approximation to the solution.

164:     Input Parameters:
165: .   user - user-defined application context
166: .   X    - vector

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

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

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


200: /* ------------------------------------------------------------------ */
201: /*
202:    FormFunctionGradient - Evaluates the function and corresponding gradient.

204:    Input Parameters:
205:    tao - the Tao context
206:    X   - the input vector
207:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

209:    Output Parameters:
210:    f   - the newly evaluated function
211:    G   - the newly evaluated gradient
212: */
213: PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr){

215:   AppCtx         *user = (AppCtx *)ptr;
217:   PetscInt       i,j,k,ind;
218:   PetscInt       xe,ye,xsm,ysm,xep,yep;
219:   PetscInt       xs, ys, xm, ym, gxm, gym, gxs, gys;
220:   PetscInt       mx = user->mx, my = user->my;
221:   PetscReal      three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three;
222:   PetscReal      p5 = 0.5, area, val, flin, fquad;
223:   PetscReal      v,vb,vl,vr,vt,dvdx,dvdy;
224:   PetscReal      hx = 1.0/(user->mx + 1);
225:   PetscReal      hy = 1.0/(user->my + 1);
226:   Vec            localX = user->localX;

229:   /* Initialize */
230:   flin = fquad = zero;

232:   VecSet(G, zero);
233:   /*
234:      Scatter ghost points to local vector,using the 2-step process
235:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
236:      By placing code between these two statements, computations can be
237:      done while messages are in transition.
238:   */
239:   DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX);
240:   DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX);

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

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

249:   /* Set local loop dimensions */
250:   xe = xs+xm;
251:   ye = ys+ym;
252:   if (xs == 0)  xsm = xs-1;
253:   else          xsm = xs;
254:   if (ys == 0)  ysm = ys-1;
255:   else          ysm = ys;
256:   if (xe == mx) xep = xe+1;
257:   else          xep = xe;
258:   if (ye == my) yep = ye+1;
259:   else          yep = ye;

261:   /* Compute local gradient contributions over the lower triangular elements */
262:   for (j=ysm; j<ye; j++) {  /*  for (j=-1; j<my; j++) */
263:     for (i=xsm; i<xe; i++) {  /*   for (i=-1; i<mx; i++) */
264:       k = (j-gys)*gxm + i-gxs;
265:       v = zero;
266:       vr = zero;
267:       vt = zero;
268:       if (i >= 0 && j >= 0) v = x[k];
269:       if (i < mx-1 && j > -1) vr = x[k+1];
270:       if (i > -1 && j < my-1) vt = x[k+gxm];
271:       dvdx = (vr-v)/hx;
272:       dvdy = (vt-v)/hy;
273:       if (i != -1 && j != -1) {
274:         ind = k; val = - dvdx/hx - dvdy/hy - cdiv3;
275:         VecSetValuesLocal(G,1,&k,&val,ADD_VALUES);
276:       }
277:       if (i != mx-1 && j != -1) {
278:         ind = k+1; val =  dvdx/hx - cdiv3;
279:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
280:       }
281:       if (i != -1 && j != my-1) {
282:         ind = k+gxm; val = dvdy/hy - cdiv3;
283:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
284:       }
285:       fquad += dvdx*dvdx + dvdy*dvdy;
286:       flin -= cdiv3 * (v + vr + vt);
287:     }
288:   }

290:   /* Compute local gradient contributions over the upper triangular elements */
291:   for (j=ys; j<yep; j++) { /*  for (j=0; j<=my; j++) */
292:     for (i=xs; i<xep; i++) {  /*   for (i=0; i<=mx; i++) */
293:       k = (j-gys)*gxm + i-gxs;
294:       vb = zero;
295:       vl = zero;
296:       v  = zero;
297:       if (i < mx && j > 0) vb = x[k-gxm];
298:       if (i > 0 && j < my) vl = x[k-1];
299:       if (i < mx && j < my) v = x[k];
300:       dvdx = (v-vl)/hx;
301:       dvdy = (v-vb)/hy;
302:       if (i != mx && j != 0) {
303:         ind = k-gxm; val = - dvdy/hy - cdiv3;
304:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
305:       }
306:       if (i != 0 && j != my) {
307:         ind = k-1; val =  - dvdx/hx - cdiv3;
308:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
309:       }
310:       if (i != mx && j != my) {
311:         ind = k; val =  dvdx/hx + dvdy/hy - cdiv3;
312:         VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES);
313:       }
314:       fquad += dvdx*dvdx + dvdy*dvdy;
315:       flin -= cdiv3 * (vb + vl + v);
316:     }
317:   }


320:   /* Restore vector */
321:   VecRestoreArray(localX,&x);

323:   /* Assemble gradient vector */
324:   VecAssemblyBegin(G);
325:   VecAssemblyEnd(G);

327:   /* Scale the gradient */
328:   area = p5*hx*hy;
329:   floc = area * (p5 * fquad + flin);
330:   VecScale(G, area);

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

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

337:   return(0);
338: }



342: PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void*ctx)
343: {
344:   AppCtx         *user= (AppCtx*) ctx;
346:   PetscInt       i,j,k;
347:   PetscInt       col[5],row;
348:   PetscInt       xs,xm,gxs,gxm,ys,ym,gys,gym;
349:   PetscReal      v[5];
350:   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;

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

354:   /*
355:      Get local grid boundaries
356:   */

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

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

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

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

368:       k=0;
369:       if (j>gys){
370:         v[k]=-2*hyhy; col[k]=row - gxm; k++;
371:       }

373:       if (i>gxs){
374:         v[k]= -2*hxhx; col[k]=row - 1; k++;
375:       }

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

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

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

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

389:     }
390:   }
391:   /*
392:      Assemble matrix, using the 2-step process:
393:      MatAssemblyBegin(), MatAssemblyEnd().
394:      By placing code between these two statements, computations can be
395:      done while messages are in transition.
396:   */
397:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
398:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
399:   /*
400:     Tell the matrix we will never add a new nonzero location to the
401:     matrix. If we do it will generate an error.
402:   */
403:   MatScale(A,area);
404:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
405:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
406:   PetscLogFlops(9*xm*ym+49*xm);
407:   return(0);
408: }