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:   PetscCall(PetscInitialize(&argc, &argv, NULL, 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:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg));
 84:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg));
 85:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg));

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

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

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

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

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

104:   /* The TAO code begins here */

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

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

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

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

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

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

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

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

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

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

141:   PetscCall(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:   PetscFunctionBegin;
163:   /* Get local mesh boundaries */
164:   PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL));
165:   PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL));

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

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

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

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

210:   PetscFunctionBegin;
211:   /* Initialize */
212:   flin = fquad = zero;

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

224:   /* Get pointer to vector data */
225:   PetscCall(VecGetArray(localX, &x));

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

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

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

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

307:   /* Restore vector */
308:   PetscCall(VecRestoreArray(localX, &x));

310:   /* Assemble gradient vector */
311:   PetscCall(VecAssemblyBegin(G));
312:   PetscCall(VecAssemblyEnd(G));

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

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

322:   PetscCall(PetscLogFlops((ye - ysm) * (xe - xsm) * 20 + (xep - xs) * (yep - ys) * 16));
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

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

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

337:   /*
338:      Get local grid boundaries
339:   */

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

345:   for (j = ys; j < ys + ym; j++) {
346:     for (i = xs; i < xs + xm; i++) {
347:       row = (j - gys) * gxm + (i - gxs);

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

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

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

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

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

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

400: /*TEST

402:    build:
403:       requires: !complex

405:    test:
406:       suffix: 1
407:       nsize: 2
408:       args: -tao_monitor_short -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4

410: TEST*/