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