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