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

 44: /*
 45:    User-defined application context - contains data needed by the
 46:    application-provided call-back routines, FormFunction() and
 47:    FormGradient().
 48: */
 49: typedef struct {
 50:   /* parameters */
 51:   PetscInt  mx, my; /* global discretization in x- and y-directions */
 52:   PetscReal param;  /* nonlinearity parameter */

 54:   /* work space */
 55:   Vec localX; /* local vectors */
 56:   DM  dm;     /* distributed array data structure */
 57: } AppCtx;

 59: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
 60: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 61: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);

 63: int main(int argc, char **argv)
 64: {
 65:   Vec       x;
 66:   Mat       H;
 67:   PetscInt  Nx, Ny;
 68:   Tao       tao;
 69:   PetscBool flg;
 70:   KSP       ksp;
 71:   PC        pc;
 72:   AppCtx    user;

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

 76:   /* Specify default dimension of the problem */
 77:   user.param = 5.0;
 78:   user.mx    = 10;
 79:   user.my    = 10;
 80:   Nx = Ny = PETSC_DECIDE;

 82:   /* Check for any command line arguments that override defaults */
 83:   PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg);
 84:   PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg);
 85:   PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg);

 87:   PetscPrintf(PETSC_COMM_WORLD, "\n---- Elastic-Plastic Torsion Problem -----\n");
 88:   PetscPrintf(PETSC_COMM_WORLD, "mx: %" PetscInt_FMT "     my: %" PetscInt_FMT "   \n\n", user.mx, user.my);

 90:   /* Set up distributed array */
 91:   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);
 92:   DMSetFromOptions(user.dm);
 93:   DMSetUp(user.dm);

 95:   /* Create vectors */
 96:   DMCreateGlobalVector(user.dm, &x);

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

100:   /* Create Hessian */
101:   DMCreateMatrix(user.dm, &H);
102:   MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE);

104:   /* The TAO code begins here */

106:   /* Create TAO solver and set desired solution method */
107:   TaoCreate(PETSC_COMM_WORLD, &tao);
108:   TaoSetType(tao, TAOCG);

110:   /* Set initial solution guess */
111:   FormInitialGuess(&user, x);
112:   TaoSetSolution(tao, x);

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

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

119:   /* Check for any TAO command line options */
120:   TaoSetFromOptions(tao);

122:   TaoGetKSP(tao, &ksp);
123:   if (ksp) {
124:     KSPGetPC(ksp, &pc);
125:     PCSetType(pc, PCNONE);
126:   }

128:   /* SOLVE THE APPLICATION */
129:   TaoSolve(tao);

131:   /* Free TAO data structures */
132:   TaoDestroy(&tao);

134:   /* Free PETSc data structures */
135:   VecDestroy(&x);
136:   MatDestroy(&H);

138:   VecDestroy(&user.localX);
139:   DMDestroy(&user.dm);

141:   PetscFinalize();
142:   return 0;
143: }

145: /* ------------------------------------------------------------------- */
146: /*
147:     FormInitialGuess - Computes an initial approximation to the solution.

149:     Input Parameters:
150: .   user - user-defined application context
151: .   X    - vector

153:     Output Parameters:
154:     X    - vector
155: */
156: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
157: {
158:   PetscInt  i, j, k, mx = user->mx, my = user->my;
159:   PetscInt  xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye;
160:   PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), temp, val;

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

166:   /* Compute initial guess over locally owned part of mesh */
167:   xe = xs + xm;
168:   ye = ys + ym;
169:   for (j = ys; j < ye; j++) { /*  for (j=0; j<my; j++) */
170:     temp = PetscMin(j + 1, my - j) * hy;
171:     for (i = xs; i < xe; i++) { /*  for (i=0; i<mx; i++) */
172:       k   = (j - gys) * gxm + i - gxs;
173:       val = PetscMin((PetscMin(i + 1, mx - i)) * hx, temp);
174:       VecSetValuesLocal(X, 1, &k, &val, ADD_VALUES);
175:     }
176:   }
177:   VecAssemblyBegin(X);
178:   VecAssemblyEnd(X);
179:   return 0;
180: }

182: /* ------------------------------------------------------------------ */
183: /*
184:    FormFunctionGradient - Evaluates the function and corresponding gradient.

186:    Input Parameters:
187:    tao - the Tao context
188:    X   - the input vector
189:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

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

209:   /* Initialize */
210:   flin = fquad = zero;

212:   VecSet(G, zero);
213:   /*
214:      Scatter ghost points to local vector,using the 2-step process
215:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
216:      By placing code between these two statements, computations can be
217:      done while messages are in transition.
218:   */
219:   DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
220:   DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);

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

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

229:   /* Set local loop dimensions */
230:   xe = xs + xm;
231:   ye = ys + ym;
232:   if (xs == 0) xsm = xs - 1;
233:   else xsm = xs;
234:   if (ys == 0) ysm = ys - 1;
235:   else ysm = ys;
236:   if (xe == mx) xep = xe + 1;
237:   else xep = xe;
238:   if (ye == my) yep = ye + 1;
239:   else yep = ye;

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

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

305:   /* Restore vector */
306:   VecRestoreArray(localX, &x);

308:   /* Assemble gradient vector */
309:   VecAssemblyBegin(G);
310:   VecAssemblyEnd(G);

312:   /* Scale the gradient */
313:   area = p5 * hx * hy;
314:   floc = area * (p5 * fquad + flin);
315:   VecScale(G, area);

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

320:   PetscLogFlops((ye - ysm) * (xe - xsm) * 20 + (xep - xs) * (yep - ys) * 16);
321:   return 0;
322: }

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

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

335:   /*
336:      Get local grid boundaries
337:   */

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

342:   for (j = ys; j < ys + ym; j++) {
343:     for (i = xs; i < xs + xm; i++) {
344:       row = (j - gys) * gxm + (i - gxs);

346:       k = 0;
347:       if (j > gys) {
348:         v[k]   = -2 * hyhy;
349:         col[k] = row - gxm;
350:         k++;
351:       }

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

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

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

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

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

397: /*TEST

399:    build:
400:       requires: !complex

402:    test:
403:       suffix: 1
404:       nsize: 2
405:       args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4

407: TEST*/