1 /* Program usage: mpiexec -n <proc> eptorsion2 [-help] [all TAO options] */ 2 3 /* ---------------------------------------------------------------------- 4 5 Elastic-plastic torsion problem. 6 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: 10 11 min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)} 12 13 where C is the torsion angle per unit length. 14 15 The uniprocessor version of this code is eptorsion1.c; the Fortran 16 version of this code is eptorsion2f.F. 17 18 This application solves the problem without calculating hessians 19 ---------------------------------------------------------------------- */ 20 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 */ 32 33 #include <petsctao.h> 34 #include <petscdmda.h> 35 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"; 43 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 */ 53 54 /* work space */ 55 Vec localX; /* local vectors */ 56 DM dm; /* distributed array data structure */ 57 } AppCtx; 58 59 PetscErrorCode FormInitialGuess(AppCtx *, Vec); 60 PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *); 61 PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *); 62 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; 73 74 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 75 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; 81 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)); 86 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)); 89 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)); 94 95 /* Create vectors */ 96 PetscCall(DMCreateGlobalVector(user.dm, &x)); 97 98 PetscCall(DMCreateLocalVector(user.dm, &user.localX)); 99 100 /* Create Hessian */ 101 PetscCall(DMCreateMatrix(user.dm, &H)); 102 PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE)); 103 104 /* The TAO code begins here */ 105 106 /* Create TAO solver and set desired solution method */ 107 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 108 PetscCall(TaoSetType(tao, TAOCG)); 109 110 /* Set initial solution guess */ 111 PetscCall(FormInitialGuess(&user, x)); 112 PetscCall(TaoSetSolution(tao, x)); 113 114 /* Set routine for function and gradient evaluation */ 115 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user)); 116 117 PetscCall(TaoSetHessian(tao, H, H, FormHessian, (void *)&user)); 118 119 /* Check for any TAO command line options */ 120 PetscCall(TaoSetFromOptions(tao)); 121 122 PetscCall(TaoGetKSP(tao, &ksp)); 123 if (ksp) { 124 PetscCall(KSPGetPC(ksp, &pc)); 125 PetscCall(PCSetType(pc, PCNONE)); 126 } 127 128 /* SOLVE THE APPLICATION */ 129 PetscCall(TaoSolve(tao)); 130 131 /* Free TAO data structures */ 132 PetscCall(TaoDestroy(&tao)); 133 134 /* Free PETSc data structures */ 135 PetscCall(VecDestroy(&x)); 136 PetscCall(MatDestroy(&H)); 137 138 PetscCall(VecDestroy(&user.localX)); 139 PetscCall(DMDestroy(&user.dm)); 140 141 PetscCall(PetscFinalize()); 142 return 0; 143 } 144 145 /* ------------------------------------------------------------------- */ 146 /* 147 FormInitialGuess - Computes an initial approximation to the solution. 148 149 Input Parameters: 150 . user - user-defined application context 151 . X - vector 152 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; 161 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)); 166 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 } 182 183 /* ------------------------------------------------------------------ */ 184 /* 185 FormFunctionGradient - Evaluates the function and corresponding gradient. 186 187 Input Parameters: 188 tao - the Tao context 189 X - the input vector 190 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 191 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; 209 210 PetscFunctionBegin; 211 /* Initialize */ 212 flin = fquad = zero; 213 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)); 223 224 /* Get pointer to vector data */ 225 PetscCall(VecGetArray(localX, &x)); 226 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)); 230 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; 242 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 } 274 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 } 306 307 /* Restore vector */ 308 PetscCall(VecRestoreArray(localX, &x)); 309 310 /* Assemble gradient vector */ 311 PetscCall(VecAssemblyBegin(G)); 312 PetscCall(VecAssemblyEnd(G)); 313 314 /* Scale the gradient */ 315 area = p5 * hx * hy; 316 floc = area * (p5 * fquad + flin); 317 PetscCall(VecScale(G, area)); 318 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)); 321 322 PetscCall(PetscLogFlops((ye - ysm) * (xe - xsm) * 20 + (xep - xs) * (yep - ys) * 16)); 323 PetscFunctionReturn(PETSC_SUCCESS); 324 } 325 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; 334 335 /* Compute the quadratic term in the objective function */ 336 337 /* 338 Get local grid boundaries 339 */ 340 341 PetscFunctionBegin; 342 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 343 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 344 345 for (j = ys; j < ys + ym; j++) { 346 for (i = xs; i < xs + xm; i++) { 347 row = (j - gys) * gxm + (i - gxs); 348 349 k = 0; 350 if (j > gys) { 351 v[k] = -2 * hyhy; 352 col[k] = row - gxm; 353 k++; 354 } 355 356 if (i > gxs) { 357 v[k] = -2 * hxhx; 358 col[k] = row - 1; 359 k++; 360 } 361 362 v[k] = 4.0 * (hxhx + hyhy); 363 col[k] = row; 364 k++; 365 366 if (i + 1 < gxs + gxm) { 367 v[k] = -2.0 * hxhx; 368 col[k] = row + 1; 369 k++; 370 } 371 372 if (j + 1 < gys + gym) { 373 v[k] = -2 * hyhy; 374 col[k] = row + gxm; 375 k++; 376 } 377 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 } 399 400 /*TEST 401 402 build: 403 requires: !complex 404 405 test: 406 suffix: 1 407 nsize: 2 408 args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4 409 410 TEST*/ 411