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 Vec x; 65 Mat H; 66 PetscInt Nx, Ny; 67 Tao tao; 68 PetscBool flg; 69 KSP ksp; 70 PC pc; 71 AppCtx user; 72 73 PetscInitialize(&argc, &argv, (char *)0, help); 74 75 /* Specify default dimension of the problem */ 76 user.param = 5.0; 77 user.mx = 10; 78 user.my = 10; 79 Nx = Ny = PETSC_DECIDE; 80 81 /* Check for any command line arguments that override defaults */ 82 PetscCall(PetscOptionsGetReal(NULL, NULL, "-par", &user.param, &flg)); 83 PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &user.my, &flg)); 84 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &user.mx, &flg)); 85 86 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n---- Elastic-Plastic Torsion Problem -----\n")); 87 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mx: %" PetscInt_FMT " my: %" PetscInt_FMT " \n\n", user.mx, user.my)); 88 89 /* Set up distributed array */ 90 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)); 91 PetscCall(DMSetFromOptions(user.dm)); 92 PetscCall(DMSetUp(user.dm)); 93 94 /* Create vectors */ 95 PetscCall(DMCreateGlobalVector(user.dm, &x)); 96 97 PetscCall(DMCreateLocalVector(user.dm, &user.localX)); 98 99 /* Create Hessian */ 100 PetscCall(DMCreateMatrix(user.dm, &H)); 101 PetscCall(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE)); 102 103 /* The TAO code begins here */ 104 105 /* Create TAO solver and set desired solution method */ 106 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao)); 107 PetscCall(TaoSetType(tao, TAOCG)); 108 109 /* Set initial solution guess */ 110 PetscCall(FormInitialGuess(&user, x)); 111 PetscCall(TaoSetSolution(tao, x)); 112 113 /* Set routine for function and gradient evaluation */ 114 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user)); 115 116 PetscCall(TaoSetHessian(tao, H, H, FormHessian, (void *)&user)); 117 118 /* Check for any TAO command line options */ 119 PetscCall(TaoSetFromOptions(tao)); 120 121 PetscCall(TaoGetKSP(tao, &ksp)); 122 if (ksp) { 123 PetscCall(KSPGetPC(ksp, &pc)); 124 PetscCall(PCSetType(pc, PCNONE)); 125 } 126 127 /* SOLVE THE APPLICATION */ 128 PetscCall(TaoSolve(tao)); 129 130 /* Free TAO data structures */ 131 PetscCall(TaoDestroy(&tao)); 132 133 /* Free PETSc data structures */ 134 PetscCall(VecDestroy(&x)); 135 PetscCall(MatDestroy(&H)); 136 137 PetscCall(VecDestroy(&user.localX)); 138 PetscCall(DMDestroy(&user.dm)); 139 140 PetscFinalize(); 141 return 0; 142 } 143 144 /* ------------------------------------------------------------------- */ 145 /* 146 FormInitialGuess - Computes an initial approximation to the solution. 147 148 Input Parameters: 149 . user - user-defined application context 150 . X - vector 151 152 Output Parameters: 153 X - vector 154 */ 155 PetscErrorCode FormInitialGuess(AppCtx *user, Vec X) { 156 PetscInt i, j, k, mx = user->mx, my = user->my; 157 PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye; 158 PetscReal hx = 1.0 / (mx + 1), hy = 1.0 / (my + 1), temp, val; 159 160 PetscFunctionBegin; 161 /* Get local mesh boundaries */ 162 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 163 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 164 165 /* Compute initial guess over locally owned part of mesh */ 166 xe = xs + xm; 167 ye = ys + ym; 168 for (j = ys; j < ye; j++) { /* for (j=0; j<my; j++) */ 169 temp = PetscMin(j + 1, my - j) * hy; 170 for (i = xs; i < xe; i++) { /* for (i=0; i<mx; i++) */ 171 k = (j - gys) * gxm + i - gxs; 172 val = PetscMin((PetscMin(i + 1, mx - i)) * hx, temp); 173 PetscCall(VecSetValuesLocal(X, 1, &k, &val, ADD_VALUES)); 174 } 175 } 176 PetscCall(VecAssemblyBegin(X)); 177 PetscCall(VecAssemblyEnd(X)); 178 PetscFunctionReturn(0); 179 } 180 181 /* ------------------------------------------------------------------ */ 182 /* 183 FormFunctionGradient - Evaluates the function and corresponding gradient. 184 185 Input Parameters: 186 tao - the Tao context 187 X - the input vector 188 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 189 190 Output Parameters: 191 f - the newly evaluated function 192 G - the newly evaluated gradient 193 */ 194 PetscErrorCode FormFunctionGradient(Tao tao, Vec X, PetscReal *f, Vec G, void *ptr) { 195 AppCtx *user = (AppCtx *)ptr; 196 PetscInt i, j, k, ind; 197 PetscInt xe, ye, xsm, ysm, xep, yep; 198 PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys; 199 PetscInt mx = user->mx, my = user->my; 200 PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param / three; 201 PetscReal p5 = 0.5, area, val, flin, fquad; 202 PetscReal v, vb, vl, vr, vt, dvdx, dvdy; 203 PetscReal hx = 1.0 / (user->mx + 1); 204 PetscReal hy = 1.0 / (user->my + 1); 205 Vec localX = user->localX; 206 207 PetscFunctionBegin; 208 /* Initialize */ 209 flin = fquad = zero; 210 211 PetscCall(VecSet(G, zero)); 212 /* 213 Scatter ghost points to local vector,using the 2-step process 214 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 215 By placing code between these two statements, computations can be 216 done while messages are in transition. 217 */ 218 PetscCall(DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX)); 219 PetscCall(DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX)); 220 221 /* Get pointer to vector data */ 222 PetscCall(VecGetArray(localX, &x)); 223 224 /* Get local mesh boundaries */ 225 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 226 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 227 228 /* Set local loop dimensions */ 229 xe = xs + xm; 230 ye = ys + ym; 231 if (xs == 0) xsm = xs - 1; 232 else xsm = xs; 233 if (ys == 0) ysm = ys - 1; 234 else ysm = ys; 235 if (xe == mx) xep = xe + 1; 236 else xep = xe; 237 if (ye == my) yep = ye + 1; 238 else yep = ye; 239 240 /* Compute local gradient contributions over the lower triangular elements */ 241 for (j = ysm; j < ye; j++) { /* for (j=-1; j<my; j++) */ 242 for (i = xsm; i < xe; i++) { /* for (i=-1; i<mx; i++) */ 243 k = (j - gys) * gxm + i - gxs; 244 v = zero; 245 vr = zero; 246 vt = zero; 247 if (i >= 0 && j >= 0) v = x[k]; 248 if (i < mx - 1 && j > -1) vr = x[k + 1]; 249 if (i > -1 && j < my - 1) vt = x[k + gxm]; 250 dvdx = (vr - v) / hx; 251 dvdy = (vt - v) / hy; 252 if (i != -1 && j != -1) { 253 ind = k; 254 val = -dvdx / hx - dvdy / hy - cdiv3; 255 PetscCall(VecSetValuesLocal(G, 1, &k, &val, ADD_VALUES)); 256 } 257 if (i != mx - 1 && j != -1) { 258 ind = k + 1; 259 val = dvdx / hx - cdiv3; 260 PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES)); 261 } 262 if (i != -1 && j != my - 1) { 263 ind = k + gxm; 264 val = dvdy / hy - cdiv3; 265 PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES)); 266 } 267 fquad += dvdx * dvdx + dvdy * dvdy; 268 flin -= cdiv3 * (v + vr + vt); 269 } 270 } 271 272 /* Compute local gradient contributions over the upper triangular elements */ 273 for (j = ys; j < yep; j++) { /* for (j=0; j<=my; j++) */ 274 for (i = xs; i < xep; i++) { /* for (i=0; i<=mx; i++) */ 275 k = (j - gys) * gxm + i - gxs; 276 vb = zero; 277 vl = zero; 278 v = zero; 279 if (i < mx && j > 0) vb = x[k - gxm]; 280 if (i > 0 && j < my) vl = x[k - 1]; 281 if (i < mx && j < my) v = x[k]; 282 dvdx = (v - vl) / hx; 283 dvdy = (v - vb) / hy; 284 if (i != mx && j != 0) { 285 ind = k - gxm; 286 val = -dvdy / hy - cdiv3; 287 PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES)); 288 } 289 if (i != 0 && j != my) { 290 ind = k - 1; 291 val = -dvdx / hx - cdiv3; 292 PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES)); 293 } 294 if (i != mx && j != my) { 295 ind = k; 296 val = dvdx / hx + dvdy / hy - cdiv3; 297 PetscCall(VecSetValuesLocal(G, 1, &ind, &val, ADD_VALUES)); 298 } 299 fquad += dvdx * dvdx + dvdy * dvdy; 300 flin -= cdiv3 * (vb + vl + v); 301 } 302 } 303 304 /* Restore vector */ 305 PetscCall(VecRestoreArray(localX, &x)); 306 307 /* Assemble gradient vector */ 308 PetscCall(VecAssemblyBegin(G)); 309 PetscCall(VecAssemblyEnd(G)); 310 311 /* Scale the gradient */ 312 area = p5 * hx * hy; 313 floc = area * (p5 * fquad + flin); 314 PetscCall(VecScale(G, area)); 315 316 /* Sum function contributions from all processes */ /* TODO: Change to PetscCallMPI() */ 317 PetscCall((PetscErrorCode)MPI_Allreduce((void *)&floc, (void *)f, 1, MPIU_REAL, MPIU_SUM, MPI_COMM_WORLD)); 318 319 PetscCall(PetscLogFlops((ye - ysm) * (xe - xsm) * 20 + (xep - xs) * (yep - ys) * 16)); 320 PetscFunctionReturn(0); 321 } 322 323 PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void *ctx) { 324 AppCtx *user = (AppCtx *)ctx; 325 PetscInt i, j, k; 326 PetscInt col[5], row; 327 PetscInt xs, xm, gxs, gxm, ys, ym, gys, gym; 328 PetscReal v[5]; 329 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; 330 331 /* Compute the quadratic term in the objective function */ 332 333 /* 334 Get local grid boundaries 335 */ 336 337 PetscFunctionBegin; 338 PetscCall(DMDAGetCorners(user->dm, &xs, &ys, NULL, &xm, &ym, NULL)); 339 PetscCall(DMDAGetGhostCorners(user->dm, &gxs, &gys, NULL, &gxm, &gym, NULL)); 340 341 for (j = ys; j < ys + ym; j++) { 342 for (i = xs; i < xs + xm; i++) { 343 row = (j - gys) * gxm + (i - gxs); 344 345 k = 0; 346 if (j > gys) { 347 v[k] = -2 * hyhy; 348 col[k] = row - gxm; 349 k++; 350 } 351 352 if (i > gxs) { 353 v[k] = -2 * hxhx; 354 col[k] = row - 1; 355 k++; 356 } 357 358 v[k] = 4.0 * (hxhx + hyhy); 359 col[k] = row; 360 k++; 361 362 if (i + 1 < gxs + gxm) { 363 v[k] = -2.0 * hxhx; 364 col[k] = row + 1; 365 k++; 366 } 367 368 if (j + 1 < gys + gym) { 369 v[k] = -2 * hyhy; 370 col[k] = row + gxm; 371 k++; 372 } 373 374 PetscCall(MatSetValuesLocal(A, 1, &row, k, col, v, INSERT_VALUES)); 375 } 376 } 377 /* 378 Assemble matrix, using the 2-step process: 379 MatAssemblyBegin(), MatAssemblyEnd(). 380 By placing code between these two statements, computations can be 381 done while messages are in transition. 382 */ 383 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 384 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 385 /* 386 Tell the matrix we will never add a new nonzero location to the 387 matrix. If we do it will generate an error. 388 */ 389 PetscCall(MatScale(A, area)); 390 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 391 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE)); 392 PetscCall(PetscLogFlops(9 * xm * ym + 49 * xm)); 393 PetscFunctionReturn(0); 394 } 395 396 /*TEST 397 398 build: 399 requires: !complex 400 401 test: 402 suffix: 1 403 nsize: 2 404 args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4 405 406 TEST*/ 407