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, PetscCtx 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*/