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