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 /* 46 User-defined application context - contains data needed by the 47 application-provided call-back routines, FormFunction() and 48 FormGradient(). 49 */ 50 typedef struct { 51 /* parameters */ 52 PetscInt mx, my; /* global discretization in x- and y-directions */ 53 PetscReal param; /* nonlinearity parameter */ 54 55 /* work space */ 56 Vec localX; /* local vectors */ 57 DM dm; /* distributed array data structure */ 58 } AppCtx; 59 60 PetscErrorCode FormInitialGuess(AppCtx*, Vec); 61 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 62 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*); 63 64 int main(int argc, char **argv) 65 { 66 Vec x; 67 Mat H; 68 PetscInt Nx, Ny; 69 Tao tao; 70 PetscBool flg; 71 KSP ksp; 72 PC pc; 73 AppCtx user; 74 75 PetscInitialize(&argc, &argv, (char *)0, help); 76 77 /* Specify default dimension of the problem */ 78 user.param = 5.0; user.mx = 10; 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: %D my: %D \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 { 157 PetscInt i, j, k, mx = user->mx, my = user->my; 158 PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys, xe, ye; 159 PetscReal hx = 1.0/(mx+1), hy = 1.0/(my+1), temp, val; 160 161 PetscFunctionBegin; 162 /* Get local mesh boundaries */ 163 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 164 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 165 166 /* Compute initial guess over locally owned part of mesh */ 167 xe = xs+xm; 168 ye = ys+ym; 169 for (j=ys; j<ye; j++) { /* for (j=0; j<my; j++) */ 170 temp = PetscMin(j+1,my-j)*hy; 171 for (i=xs; i<xe; i++) { /* for (i=0; i<mx; i++) */ 172 k = (j-gys)*gxm + i-gxs; 173 val = PetscMin((PetscMin(i+1,mx-i))*hx,temp); 174 PetscCall(VecSetValuesLocal(X,1,&k,&val,ADD_VALUES)); 175 } 176 } 177 PetscCall(VecAssemblyBegin(X)); 178 PetscCall(VecAssemblyEnd(X)); 179 PetscFunctionReturn(0); 180 } 181 182 /* ------------------------------------------------------------------ */ 183 /* 184 FormFunctionGradient - Evaluates the function and corresponding gradient. 185 186 Input Parameters: 187 tao - the Tao context 188 X - the input vector 189 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 190 191 Output Parameters: 192 f - the newly evaluated function 193 G - the newly evaluated gradient 194 */ 195 PetscErrorCode FormFunctionGradient(Tao tao,Vec X,PetscReal *f,Vec G,void *ptr) 196 { 197 AppCtx *user = (AppCtx *)ptr; 198 PetscInt i,j,k,ind; 199 PetscInt xe,ye,xsm,ysm,xep,yep; 200 PetscInt xs, ys, xm, ym, gxm, gym, gxs, gys; 201 PetscInt mx = user->mx, my = user->my; 202 PetscReal three = 3.0, zero = 0.0, *x, floc, cdiv3 = user->param/three; 203 PetscReal p5 = 0.5, area, val, flin, fquad; 204 PetscReal v,vb,vl,vr,vt,dvdx,dvdy; 205 PetscReal hx = 1.0/(user->mx + 1); 206 PetscReal hy = 1.0/(user->my + 1); 207 Vec localX = user->localX; 208 209 PetscFunctionBegin; 210 /* Initialize */ 211 flin = fquad = zero; 212 213 PetscCall(VecSet(G, zero)); 214 /* 215 Scatter ghost points to local vector,using the 2-step process 216 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 217 By placing code between these two statements, computations can be 218 done while messages are in transition. 219 */ 220 PetscCall(DMGlobalToLocalBegin(user->dm,X,INSERT_VALUES,localX)); 221 PetscCall(DMGlobalToLocalEnd(user->dm,X,INSERT_VALUES,localX)); 222 223 /* Get pointer to vector data */ 224 PetscCall(VecGetArray(localX,&x)); 225 226 /* Get local mesh boundaries */ 227 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 228 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 229 230 /* Set local loop dimensions */ 231 xe = xs+xm; 232 ye = ys+ym; 233 if (xs == 0) xsm = xs-1; 234 else xsm = xs; 235 if (ys == 0) ysm = ys-1; 236 else ysm = ys; 237 if (xe == mx) xep = xe+1; 238 else xep = xe; 239 if (ye == my) yep = ye+1; 240 else yep = ye; 241 242 /* Compute local gradient contributions over the lower triangular elements */ 243 for (j=ysm; j<ye; j++) { /* for (j=-1; j<my; j++) */ 244 for (i=xsm; i<xe; i++) { /* for (i=-1; i<mx; i++) */ 245 k = (j-gys)*gxm + i-gxs; 246 v = zero; 247 vr = zero; 248 vt = zero; 249 if (i >= 0 && j >= 0) v = x[k]; 250 if (i < mx-1 && j > -1) vr = x[k+1]; 251 if (i > -1 && j < my-1) vt = x[k+gxm]; 252 dvdx = (vr-v)/hx; 253 dvdy = (vt-v)/hy; 254 if (i != -1 && j != -1) { 255 ind = k; val = - dvdx/hx - dvdy/hy - cdiv3; 256 PetscCall(VecSetValuesLocal(G,1,&k,&val,ADD_VALUES)); 257 } 258 if (i != mx-1 && j != -1) { 259 ind = k+1; 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; val = dvdy/hy - cdiv3; 264 PetscCall(VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES)); 265 } 266 fquad += dvdx*dvdx + dvdy*dvdy; 267 flin -= cdiv3 * (v + vr + vt); 268 } 269 } 270 271 /* Compute local gradient contributions over the upper triangular elements */ 272 for (j=ys; j<yep; j++) { /* for (j=0; j<=my; j++) */ 273 for (i=xs; i<xep; i++) { /* for (i=0; i<=mx; i++) */ 274 k = (j-gys)*gxm + i-gxs; 275 vb = zero; 276 vl = zero; 277 v = zero; 278 if (i < mx && j > 0) vb = x[k-gxm]; 279 if (i > 0 && j < my) vl = x[k-1]; 280 if (i < mx && j < my) v = x[k]; 281 dvdx = (v-vl)/hx; 282 dvdy = (v-vb)/hy; 283 if (i != mx && j != 0) { 284 ind = k-gxm; val = - dvdy/hy - cdiv3; 285 PetscCall(VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES)); 286 } 287 if (i != 0 && j != my) { 288 ind = k-1; val = - dvdx/hx - cdiv3; 289 PetscCall(VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES)); 290 } 291 if (i != mx && j != my) { 292 ind = k; val = dvdx/hx + dvdy/hy - cdiv3; 293 PetscCall(VecSetValuesLocal(G,1,&ind,&val,ADD_VALUES)); 294 } 295 fquad += dvdx*dvdx + dvdy*dvdy; 296 flin -= cdiv3 * (vb + vl + v); 297 } 298 } 299 300 /* Restore vector */ 301 PetscCall(VecRestoreArray(localX,&x)); 302 303 /* Assemble gradient vector */ 304 PetscCall(VecAssemblyBegin(G)); 305 PetscCall(VecAssemblyEnd(G)); 306 307 /* Scale the gradient */ 308 area = p5*hx*hy; 309 floc = area * (p5 * fquad + flin); 310 PetscCall(VecScale(G, area)); 311 312 /* Sum function contributions from all processes */ /* TODO: Change to PetscCallMPI() */ 313 PetscCall((PetscErrorCode)MPI_Allreduce((void*)&floc,(void*)f,1,MPIU_REAL,MPIU_SUM,MPI_COMM_WORLD)); 314 315 PetscCall(PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16)); 316 PetscFunctionReturn(0); 317 } 318 319 PetscErrorCode FormHessian(Tao tao, Vec X, Mat A, Mat Hpre, void*ctx) 320 { 321 AppCtx *user= (AppCtx*) ctx; 322 PetscInt i,j,k; 323 PetscInt col[5],row; 324 PetscInt xs,xm,gxs,gxm,ys,ym,gys,gym; 325 PetscReal v[5]; 326 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; 327 328 /* Compute the quadratic term in the objective function */ 329 330 /* 331 Get local grid boundaries 332 */ 333 334 PetscFunctionBegin; 335 PetscCall(DMDAGetCorners(user->dm,&xs,&ys,NULL,&xm,&ym,NULL)); 336 PetscCall(DMDAGetGhostCorners(user->dm,&gxs,&gys,NULL,&gxm,&gym,NULL)); 337 338 for (j=ys; j<ys+ym; j++) { 339 340 for (i=xs; i< xs+xm; i++) { 341 342 row=(j-gys)*gxm + (i-gxs); 343 344 k=0; 345 if (j>gys) { 346 v[k]=-2*hyhy; col[k]=row - gxm; k++; 347 } 348 349 if (i>gxs) { 350 v[k]= -2*hxhx; col[k]=row - 1; k++; 351 } 352 353 v[k]= 4.0*(hxhx+hyhy); col[k]=row; k++; 354 355 if (i+1 < gxs+gxm) { 356 v[k]= -2.0*hxhx; col[k]=row+1; k++; 357 } 358 359 if (j+1 <gys+gym) { 360 v[k]= -2*hyhy; col[k] = row+gxm; k++; 361 } 362 363 PetscCall(MatSetValuesLocal(A,1,&row,k,col,v,INSERT_VALUES)); 364 365 } 366 } 367 /* 368 Assemble matrix, using the 2-step process: 369 MatAssemblyBegin(), MatAssemblyEnd(). 370 By placing code between these two statements, computations can be 371 done while messages are in transition. 372 */ 373 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 374 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 375 /* 376 Tell the matrix we will never add a new nonzero location to the 377 matrix. If we do it will generate an error. 378 */ 379 PetscCall(MatScale(A,area)); 380 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 381 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 382 PetscCall(PetscLogFlops(9*xm*ym+49*xm)); 383 PetscFunctionReturn(0); 384 } 385 386 /*TEST 387 388 build: 389 requires: !complex 390 391 test: 392 suffix: 1 393 nsize: 2 394 args: -tao_smonitor -tao_type nls -mx 16 -my 16 -tao_gatol 1.e-4 395 396 TEST*/ 397